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ABSTRACT 



A numerical technique is formulated, in a computer program U2DIIF, for the 
solution of flow over an airfoil executing an arbitrary unsteady motion in an inviscid 
and incompressible medium. The technique extends the well known Panel Methods for 
steady flow into solving a non-linear unsteady flow problem arising from the 
continuous vortex shedding into the trailing wake due to the unsteady motion of the 
airfoil. Numerous case-runs are presented to verify U2DIIF computer code against 
other theoretical and/or numerical methods as well as in cases where limited 
experimental data are obtainable in literatures. These case-runs include airfoils 
undergoing a step change or a modified ramp change of angle-of-attack, airfoils 
executing harmonic oscillation in pitching and plunging motions and airfoils 
penetrating a sharp edge gust. 
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THESIS DISCLAIMER 

The reader is cautioned that computer programs developed in this research may 
not have been exercised for all cases of interest. While every effort has been made, 
within the time available, to ensure that the programs are free of computational and 
logic errors, they cannot be considered validated. Any application of these programs 
Avithout additional verification is at the risk of the user. 
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I. INTRODUCTION 



A. GENER.AL 

In this thesis, a numerical method is formulated and coded in a FORT FLAN 
computer program, codename U2DIIF (Unsteady 2-Dimensional Inviscid 
Incompressible Flow), to solve for the flow over an airfoil which is executing an 
unsteady time-dependent motion in an inviscid, incompressible medium. 

B. APPROACH 

The basic approach to this problem is the extension of a ver\' general and 
powerful technique, called Panel Methods, developed by Hess & Smith [Ref 1] for 
steady potential flow problems, to include the unsteady motion of the airfoil that is 
continuously shedding vorticity into the trailing wake. This vortex shedding process 
creates the non-linearity effects of the problem in that the wake vortices infuence the 
flow over the airfoil which in turn alters the vonex shedding as the airfoil proceeds in 
time. It is this very non-linearity of unsteady flow that distinguishes itself from the well 
known steady Panel Methods solution where the mathematical formulation of the 
problem results in a set of N linear equations in N unknowns which are solved easily 
with the standard Gaussian elimination algorithm. 

The unsteady flow problem is, however, deprived of this relatively easy solution 
technique. Instead, an iterative type of solution is needed for this non-linear problem. 
The correct mathematical model must therefore be formulated to describe the vortex 
shedding process that provides the mechanism for the iteration to proceed towards a 
converged set of solution in each time step. 

It is the objective of this thesis to develop a numerical computer program that 
performs this non-linear potential flow calculation which proceeds step by step in time. 
.At each time step, a complete set of potential flow solutions, inclusive of the airfoil 
pressure distribution, force and moment coeiFicients, and the trailing vortex wake 
pattern (strengths and positions of shed vortices), is obtained. 

C SCOPE 

The Panel Method of Hess & Smith, which utilises both the distributed sources 
and vorticities as panel singularities, for steady flow solution is described in Chapter II. 
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Chapter III formulates the mathematical model for the unsteady flow problem 
and its solution procedures, highlighting the essential features in solving the non-linear 
problem of unsteady flow. 

Chapter IV describes the computer program U2DIIF, its essential capabilities, 
limitations and the necessary input set-up for typical case-runs. 

The results of some of the case-runs are presented in Chapter V. They are 
compared with other theoretical and/or numerical methods as well as in cases where 
limited experimental data are obtainable in the literature. These case-runs include 
airfoils undergoing a step change or a modified ramp change in angle^of-attack, airfoils 
executing harmonic oscillation in both pitching and plunging motions and airfoils 
penetrating a sharp edge gust. 

In the concluding remarks of Chapter VI, the future development and application 
potential of this numerical method to other studies of unsteady 2-dimensional inviscid 
incompressible flow are mentioned. 
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II, STEADY FLOW PROBLEM FORMULATION 

A. FRAME OF REFERENCE 

Consider a 2-dimensional airfoil in motion with constant linear velocity “ Vqq as 
shown in Figure 2.1. Using an (x,y) coordinate system fixed on the airfoil, where the x- 
axis coincides with the chord line originating from the leading edge towards the trailing 
edge of the airfoil, the flow in this frame of reference is steady. That is to say, the fluid 
velocity and pressure in the flow field depend only on the spatial coordinates (x,y) and 
not on time. The airfoil then appears to be submerged in an onset flow whose velocity 
is V^o and making an angle of attack, a, with the x-axis (see Figure 2.1). 

B. STEADY FLOW PANEL METHODS 

1. Definition of Nodes and Panels 

The airfoil surface is divided into (n) straight-line segments, called panels, by 
(n+ 1) arbitrary chosen points, called nodes, distributed over the airfoil contour as 
shown in Figure 2.2. The panel numbering sequence starts with panel 1 on the lower 
surface at the airfoil trailing edge and proceeds clockwise around the airfoil contour so 
that the last panel (panel n) ends on the upper surface, also at the airfoil trailing edge. 

Notice that this numbering sequence dictates that the airfoil body always lies 
on the right hand side of the i^ panel as one proceeds from the i^ node to the (i+ 1)^ 
node. Also the and the (n+1)^ nodes coincide at the trailing edge. It therefore 
facilitates, as shown in Figure 2.2, the common definition of unit normal vector Hj and 
the unit tangential vector t for all panels, i.e., is directed outward from the body into 
the flow and t is directed from the i^ node to the (i+ 1)^ node. 

2. Distribution of Singularities 

Figure 2.2 also indicates that a uniform source distribution qj and a uniform 
vorticity distribution 7 are placed on the panel. The source strength q; varies from 

i 

panel to panei whereas the vorticity strength y remains the same for ail panels. This 
particular choice of singularity distributions is one of the many types of singularity 
combinations (it happened to be the pioneering one though) ever used in a wide variety 
of the so called Panel Methods. The success of representing the flow past an arbitrary 
shaped airfoil by surface singularity distributions lies in the fact that these singularity 
distributions automatically satisfy Laplace's equation, the governing flow equation for 
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Figure 2.1 Frame of Reference for Steady Flow, 
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Figure 2.2 Panel .VIethods Representation for Steady Flow. 
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inviscid incompressible flow, and the boundary condition at the far field (^o). In 
addition, the superposition principle applies to any linear homogeneous second order 
partial differential equation such as Laplace's equation. Therefore one can build up an 
overall complicated flow field by the combination of simple flows if the appropriate 
boundary conditions on the airfoil can be satisfied accurately. For our case the overall 
flow field (represented by the velocity potential <I>) can be built up by three simple 
flows, 



<D = (Poo + (pj + <Pv (eqn2.1) 

where (Poo is the potential of the onset flow, 

<Poo = ^00 (x cosa + y sina) (eqn 2.2) 

<p^ is the velocity potential of a source distribution of strength q(s) per unit length. 




In r ds 



(eqn 2.3) 



(py is the velocity potential of a vorticity distribution of strength y(s) per unit length. 




Y(s) 

27t 



6 ds 



(eqn 2.4) 



The integrals in Equations 2.3 and 2.4 are performed along the surface 
contour s and (r,6) are polar coordinates of any field point (x,y) measured from the 
airfoil surface at an arbitrary point as shown in Figure 2.3. The difficult task of 
evaluating these integrals has been greatly simplified by our singularity distributions 
postulated to represent the flow over the airfoil: that is. instead of integrating over the 
entire airfoil contour, we integrate on each panel along a straight line where q^ and y 
are constant, then sum up the effects of all panels. Equation 2.1 therefore becomes, 

= Vqo (x cosa + y sina) 4- ^ Jpanel j t ^ ® 1 
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r 7(s) 

<t> = Von (x cosa + V sina) + j — In r ds - | 0 ds 

2k 2k 



Figure 2.3 Potential Evaluation at a Field Point. 
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It can be seen from Equation 2.5 that <I> is completely defmed if the (n+ 1) 
unknowns, (j= l,2,....,n) and 7, can be calculated using a numerical technique yet to 
be described. Once the potential O is solved, the velocity can be evaluated by taking 
grad O. At thus point we introduce a definition of disturbance potential, <p, as the sum 
of potential due to both the source and vorticity distribution, 

9 = 9s + <Pv 

Equation 2.1 therefore reads, 

<D = (poo + (p (eqn 2.7) 



The total velocity vector is thus. 



V, 



total 



VO 



= V<Poo + 7(p 

= Voo + 7(p 



The pressure can be obtained from Bernoulli's Equation, 



(eqn 2.8) 




^ ^00 = 1 _ / ^total <2 

S'^PVoo' Voo 



(eqn 2.9) 



Notice that Figure 2.3 indicates that the field point lies off the airfoil surface, 
how'ever, we are interested in field points that are on the airfoil surface. In the case of 
steady flow, the expressions for and Cp are the same for field points lying on or 
off the airfoil surface. It is nevertheless not the same m unsteady flow, as will be seen 
in Chapter III, in that must include the rigid body motion of the airfoil when one 
evaluates field points on the airfoil surface. 

3. Boundary Conditions 

The boundary conditions to be satisfied include the flow tangency conditions 
and the Kutta Condition. The flow tangency conditions are satisfied at the exterior mid 
points, called control points, of all panels by taking the resultant velocity at each 
control point to have only (V^)j but. 
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(V"). = 0, i=l,2 ,n 



(eqn 2.10) 



f- 



where (V^). and (V"). are the tangential and normal components of the total velocity at 
the control point of the i^ panel due to the free stream and the velocities induced by 
the source and vorticity distributions on all the panels, j (j= 1,2,.... ,n). 

The Kutta condition postulates that the pressures on the upper and lower 
panels at the trailing edge be equal in order that the flow leaves the trailing edge 
smoothly. By using Bernoulli's equation for steady potential flow, this pressure 
equilibrium condition implies that the tangential velocities in the downstream direction 
at the l^*' and the n* panel control points must be equal. This fact is certainly 
consistent with the knowledge that when steady flow is established, the total circulation 
over the airfoil does not change if the tangential velocities are the same at the trailing 
edge panels. 

(V‘)j = -(V% (eqn 2.11) 

If one could explicitly express Equations 2.10 and 2.11 in terms of the unknowns 
qj (i= l,2,....,n) and y. the task is then reduced to solving a linear system of (n+1) 
simultaneous equations for the (n+ 1) unknowns. 

C. INFLUENCE COEFFICIENTS 

1. The Concept of Influence Coefficients 

The numerical technique employed in Panel Methods to manipulate 
equations 2.10 and 2.1 1 into an algebraic system of linear simultaneous equations 
involves the important concept of influence coefficients. An influence coefficient is 
defined as the velocity induced at a field point by a unit strength singularity (be it a 
point singularity or a distributed singularity) placed anywhere within the flow field. In 
this case, it is the unit strength singularity distribution on one panel. Recall that 
equations 2.10 and 2.11 simply require the computation of the normal and tangential 
velocity components at all the panel control points. The normal components of 
velocities are essential in satisfying flow tangency conditions while the tangential 
components of velocities are necessary for satisfying the Kutta condition as well as 
computing the pressure distribution. The procedure is thus to compute, at the i^^ panel 
control point, the velocity components induced by the source and vorticity 
distributions on all the panels, j (j= 1,2,.. ..,n), including the i^ panel itself Summation 
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of all the induced velocities, separately for the normal and tangential components, 
together with the free stream velocity components produces all the required (V'^). and 
(V%i=l,2,....,n. 

2. Notation for Influence Coefficient 

We shall adopt a consistent set of notation for the influence coefficients used 
throughout this documentation. It is so designated to permit easy recognition in that 
each influence coefficient contains all the associated information one needs. An 
influence coefficient is denoted with a superscript and two subsript as follows: 

^ pq 

where x denotes the type of singularity involved, we shall arbitrarily use A, B and C for 
the uniformly distributed source, uniformly distributed vorticity and point vortex 
respectively. The superscript s is an indicator telling which component the induced 
velocity is. The first subscript p identifies the field point where the induced velocity is 
evaluated. The second subscript q denotes the particular singularity contributing to the 
induced velocity. 

We thus define, for the steady flow problem, the following influence 
coefficients : 

• AAj : normal velocity component induced at the i* panel control point by unit 
strength source distribution on the panel. 

• A^.j : tangential velocity component induced at the i^ panel control point by 
unit strength source distribution on the panel. 

• : normal velocity component induced at the i^ panel control point by unit 
strength vorticity distribution on the panel. 

• B^j| : tangential velocity component induced at the i^ panel control point by 
unit strength vorticity distribution on the panel. 

3. Computation of Influence Coefficients 

The influence coefficients turn out to be related, not surprisingly, to the 
geometry of the airfoil and the manner in which the panels are formed. Specifically, as 
derived in iRef 2], the .A's and B's influence coefficients,^ due to uniformly distributed 
source or vorticity are functions of ; 

• The natural logarithm of the ratio of distance from the i^ panel control point 
(the field point) to the (j + 1)^ and nodes of the panel where singularities 
are distributed. 



I 



C's coefficients will be needed only for unsteady flow. 
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• The angle, in radian, subtended at the i^"^ panel control point (the field point) by 
the (j -r 1)'^ and nodes of the panel where singularities are distributed. 

• The trigonometry* angles of the i^ and panels. 

Referring to the geometrical quantities indicated in Figure 2.4, the 
expressions^ for these influence coefficients are : 



A"".. = — [ sin(0. - e.) In 

‘J 2tc ^ ^ ‘ rjj 

1 



+ cos(9j-ep Pjj ] , 



for i 32 j 



for i = j (eqn 2. 12) 




= ( sin(e. - ep pjj - cos(0. - 0j) In 

= 0 



V. + i 



r.. 

i] 






for i 3£ j 



for i = j (eqn 2.13) 



B". = ccs(8. — 0.) In ■ ^ 

2k ' ‘ J r.. 

1] 

= 0 




for i = j 



for i = j (eqn 2.14) 




^ ( cos(0. - 0j) p.j 4- sin(0j - 0p In Al±L ] ^ 
-7t r* * 



for i j 



for i = j (eqn 2.15) 



where : 



r|j + , = V ((xitij - X. + + (ymi-yj + ,)-| 
I'll = V [(xntij-Xj)^ + (ym|-yj)-] 
xm; = 

ym. = i'2(y. + y.^,) 



^Actual computation uses A'.. = 



time. 



— and B^jj = A"jj to reduce computing 
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Figure 2.4 Influence Coefllcienis due to Uniformly Distributed Singularities. 
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Y; ^ 1 - 



9. = arcian (— ^ — —) 

I V — V 






V. J_ , — V. 

6. = arctan( ' ’ ^ 



p.. = arctan( ^^^ - arctan(— 



xmj-Xj + ^ 



xnij - X. 



D. NUMERICAL SOLUTION SCHEME 
1. Rewriting the Boundary Conditions 

Using the concept of influence coefficients, the flow tangency conditions of 
Equation 2.10 can be expressed as. 



n n 

I(A";.qj| - rl B"i| + sin(a-6j) = 0 , i=l,2, n (eqn2.16) 

= I j = I 

The Kutta condition of Equation 2.11, in terms of influence coefficients, looks 



- E I A‘,i q, 1 - Y E B‘ii - Vco cos(a - 6,) 

j = i i=l 

= E t A‘„i q| 1 + V E B‘„i + Voo cos(a-e„) (eqn 2.17) 

i=l j=l 

The negative signs appearing on the left-hand-side of Equation 2.17 are a 
direct consequence of our definition of unit tangential vector. In other words, the 
tangential velocities on the lower surface panels downstream of the front stagnation 
pome have negative values. This Teaturc in fact allows one to predict the front 
stagnation point by interpolating the velocity distribution around the leading edge. 

2. Solving the Strengths of Source and Vorticity Distributions 

It is not difficult at this stage to see that if we collect the like terms in 
Equation 2.17 and expand Equation 2.16 for all i's (i= l,2,....,n), these equations 
constitute none other than a linear algebraic system of (n+ 1) equations as shown in 
the matrix Equation 2.18. 
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(eqn 2.18) 



^1,1 ^1,2 ^1,3 




'^r 






^2,1 ^2,2 ^2,3 ^2,n-M 




^2 




b. 


^3,1 ^3,2 ^3,3 •••••• ^3,n+l 




^3 


= 




^n,l ^n,2 ^n,n+l 








b„ 


^n+1,1 ' * * ' n + l,n+l 




Y 




^n+l , 



Equation 2.18 is a set of linearly independent equations which can be easily 
solved by any standard linear system solver, one of which is the well known method of 
Gaussian Elimination with Partial Pivoting. 

3. Computation of Velocity and Pressure Distribution 

Once the q^ (j= l,2,....,n) and y are solved, the velocities at all the panel 
control points can be evaluated. Only the tangential components exist since the normal 
components are already set to zeroes by the flow tangency conditions. 



V 

V 



total 

00 




i= 1,2,.. ..,n 



(eqn 2.19) 



where : 



n n 

(V‘)i = S 1 Atj qj ] + Y I B'jj + Voo cos(« - 0;) , i= 1,2 n (eqn 2.20) 

i=l j=l 

Substituting Equation 2.20 into the Cp equation (Equation 2.9), the pressure 
coefHcients at the i* panel control point is : 

(Cp)j = 1 - (V^).2 , i=l,2, n (eqn 2.21) 

4. Computation of Forces and Moments 

The two dimensional aerodynamic coefficients of lift (Cg), drag (C^) and 
pitching moment (C^^) about the leading edge are computed by integration of the 
pressure distribution assuming constant Cp exists in each panel. The computation is 
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first done bv intearatins forces in the airfoil-fixed coordinate svstem followed bv a 

•I* W W' ✓ 

rotation to the respective lift and drag directions along and perpendicular to the free 



stream (V^q) as follows ; 




i= 1 


(eqn 2.22) 


Cx = ■ I(Cp>i(yi+i-yi) 

i= 1 
n 


(eqn 2.23) 



= I (c )i 1(X; + 1 - Xj) xnij + (Yi + 1 - Y;) Y^i 1 (eqn 2.24) 

i= 1 



Cj = cosa + Cy sina 


(eqn 2.25) 

• 


Cp = C, cosa — sina 

V y X 


(eqn 2.26) 
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m. UNSTEADY FLOW PROBLEM FORMULATION 



A. OVERVIEW OF UNSTEADY FLOW MODELING 

1. Some PrevievYS 

Having fully understood the Panel Methods formulation and solution for the 
steady flow problem, one could then venture into the interesting and complicated 
unsteady flow case. In this Chapter, we shall see how we could build the time- 
dependency into the Panel Methods solution which has been proven to be an useful 
and accurate tool for steady flow. The approach in the unsteady flow problem 
formulation will proceed, in general, in a manner similar to Chapter II. However, as we 
go along, we will pick up the highlights of the essential differences (also similarity) 
between the two problems. Additional flow modeling of the vortex shedding process 
that greatly influences the numerical solution technique^ will be discussed in details. 

2. Specific Unsteady flow iModel 

Recall that in steady flow, the problem is considered solved as soon as the 
airfoil surface singularity distributions of source and vorticity (j= l,2,....,n) and y are 
determined. These (n+ 1) unknowns are, however, time dependent in unsteady flow. 
We therefore introduce a subscript k as the time-step counter; that is, we postulate to 
solve the unsteady flow problem at successive intervals of time, starting from tg = 0. At 
each time-step tj^ (k= l,2,....,oo), we represent the airfoil by surface singularity 
distributions consisting of source distribution (q.)j^ (j = l,2,....,n) and vorticity 
distribution Yj^. Again the source strengths vary from panel to panel but the vorticity 
strength remains the same for all panels. 

The overall circulation at time-step tj^ is simply Yj^ multiplied by the airfoil 
perimeter, €. Since the total circulation in a potential flow field must be preserved 
according to the Helmholtz s theorem of continuity of vorticity, any changes in the 
circulation on the airfoil surface must be manifested by an equal and opposite change 
in vorticit 3 / in the wake. We call this the vortex shedding process and postulate, as 
shown in Figure 3.1, that this shed vorticity takes place through a small straight line 
wake element attached as an additional panel to the trailing edge with uniform vorticity 
distribution We shall from now on refer to this panel as the shed vorticity panel, 

The shed vorticity panel will be established if its length A|^ and inclination 0j^, to x- 

^ Referring to the switch from a direct scheme to an iterative scheme. 
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Vortex Shedding at Time Step t,^ 



Helmholtz's theorem 



K (Yw>k + = r,^., 






Figure 3.1 Extension of Panel .Methods Representation for Unsteady Flow. 
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axis of the airfoil-fixed coordinate system, satisfy the Helmholtz's theorem as follows, 



\ (eqn3.1) 

\ “ ^k-l ~ “ ^*^'i'k-l ~ "''k^ (eqn J.2) 

where and Yj^.^ are respectively the total circulation and vorticity strengths which 
are already determined at a time-step tj^_^ before tj^. 

Let us project one time step ahead to i> we allow the shed vonicity panel 
to be detached from the trailing edge and get convected downstream as a concentrated 
free vortex, with circulation Aj^ (Yvv\ ^k-l *" ^k’ according to the resultant local 
velocity occurred at the center of the vortex core. At the same time a brand new shed 
vorticity panel is formed for the new time step and the process is repeated. Therefore 
the shed vorticity panel model provides exactly the desired communication mechanism 
to carry the solution from one time-step to another. 

We now return back to the time-step tj^ and immediately realise that as a 
result of this continuous vortex shedding, there has been a series of shedding processes 
occurred prior to tj^ that cummulated in a string of concentrated core vomces of 

strengths (r ^_2 ” (J'k-S “ ^k-2)> (^k-4“^k-s)’ forming the 

wake pattern behind the airfoil as shown in Figure 3.1 

The presence of the shed vorticity panel and the downstream resultant wake 
core vortices do influence the upstream flow in inviscid incompressible flow. In 
particular the shed vorticity panel itself depends on Yj^ to determine its distributed 
vorticity this in turn causes changes to (q.)j^ and Yj^- Moreover, the downstream 

core vortices that constitute the wake are convected under the influence of the free 
stream and the singularity distributions on the airfoil surface panels including the shed 
vorticity panel. The problem is thus seen to be coupled from this analytical 
standpoint. Putting this in simple mathematical terms, the algebriac system of 
equations (Equation 2.18), representing the flow tangency conditions and Kutta 
condition for steadv flow, are no longer linear because the coeiTicients a., in the left- 

i] 

hand-side matrix are not constants anymore. They are function of q^ and y instead. 
The presence of non-linearity is indeed what drives the solution scheme into an 
iterative type for unsteady flow. 
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3. Boundary Conditions 

We next investigate whether our unsteady flow model is sufTiciently 
represented, before we could proceed to search for a numerical iterative solution, by 
matching the unknowns with the available boundary conditions at time-step Recall 
that we have introduced three more unknowns \ and 0j^ in addition to 

(qj)j^ (j= 1,2,... .,n) and Yj^. We have, however, so far only identified an extra boundary 
condition, namely the Helmholtz's theorem (Equation 3.2) in conjunction with the flow 
tangency conditions at the n panel control points and the Kutta condition of pres'Sure 
equilibrium at the trailing edge panels. Clearly we are in deficit of two additional 
conditions before attempting further endeavour to solve the entire system. Basu and 
Hancock [Ref 3] suggested the following assumptions : 

• The shed vorticity panel is oriented in the direction of the local resultant 
velocity at the panel mid point. 

• The length of the shed vorticity panel is proportional to the magnitude of the 
resultant velocity at the panel mid point and the step size of the time-step. 

Following these assumiptions thus permits us to formulate two more boundary 
conditions as follows, 

(V X 

tanO. = — (eqn 3.3) 

^ (Uw)k 

\ = 'k-l) V I + (V„),^ 1 (eqn 3.4) 

where (U^y)|^ and (V^^)j^ are the total velocity components in x and y directions of the 
airfoil-fixed coordinate system. 

The flow tangency conditions are still, 

((V")4 = 0, i=l,2 n (eqn 3.5) 

However, the Kutta condition must now include the rates of change of 
potential at the trailing edge panels (unsteady Bernoulli's equation) which can be 
related directly to the rate of change of total circulation. By using a backward finite 
dificrencc approximation for this rate of change of total circulation, we express the 
Kutta condition as shown in Equation 3.6. 
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(eqn 3.6) 



I(V')ii,2 - [(V')J^2 = 2 [ |. (,p„ - <pj) ], = 2 ( ^ 

= ^pJkllia. 
tk-H;-! 

B. RIGID BODY MOTION AND FRAME OF REFERENCE 

Consider a rigid airfoil executing a time-dependent motion, comprising linear 
translation and angular rotation about the leading edge in an inviscid incompressible 
medium. We can describe this arbitrar>' motion at any time instant tj^ as the vector sum 
of a mean velocity ~Vqq, a time dependent translational velocity — [U(t)i + V(t) j] 
and a rotational velocity ~12(t) where i & j are unit vectors in the airfoil-fixed 
coordinate system as shown in Figure 3.2. 

If we continue, as in steady flow, to determine the flow with reference to the (x,y) 
coordinate system fixed on the airfoil, an observer sitting on this frame of reference 
sees an unsteadv stream velocity, V , made up by the vector sum of a mean 
velocity V^, a time dependent translational velocity [U(t) i ^ V(t) j] and a rotational 
velocity F2(t). Therefore in this frame of reference, unlike the previous case where the 
airfoil is allowed to move only with constant linear velocity, the flow is still unsteady in 
that is time dependent. 

^stream = ^oo + [(U(t) i + V(t) j)] + ^(t) (y i - X j) (eqn 3.7) 

We redefine our disturbance potential to include the potential contributions <p^ 
and from the shed vorticity panel and the wake core vortices respectively. Thus : 

9 = <Ps + <Pv + <Pw + ^cv 



We then write the total velocity with respect to this frame of reference as, 



V 



total 



= V 



stream 



7(p 



(eqn 3.9) 



Notice that this total velocity is obviously NOT the absolute velocity with respect 
to an inertial coordinate system. Such an inertial coordinate system will be the one 
where an observer only sees an on-coming flow of V^q with constant magnitude and 
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V 



^u•eam 



= 



~o 



- [(UlC) i - V7t) i)| - CIU) ly i - X j) 



Figure 3.2 Frame of Reference for Lnstcady Flow. 
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direction. We have to make this distinction clear because in our model on convection 
of core vortices, we break up the caiulation into two steps; we first convecc the core 
vortices using the resultant absolute velocity with respect to an inertial coordinate 
system, followed by determining their positions with coordinates relative to the airfoil- 
fixed axes. 

The unsteady flow Bernoulli's equation for the pressure coefTicients on the airfoil 
surface must be written with respect to the airfoil-fixed coordinate system also. Giesing 
[Ref 4] showed this to be written, in our notation, as : 




^ ^00 _ / ^stream \2 _ / ^total \2 _ 

'/^pVoo^ Voo ^co 



2 5(p 



(eqn 3.10) 



where V^^j-eam’ "^total defined according to Equations 3.7 and 3.9. 

Equations 3.8, 3.9, and 3.10 can be correlated to their counter-parts in steady 
flow, namely 2.6, 2.8 and 2.9 respectively with Vg^ream Equation 3.7 replacing the 
in Equation 2.8. 

C. TIME-DEPENDENT INFLUENCE COEFFICIENTS 
1. Definition of Time-Dependent Influence Coefficients 

The influence coefficients. A”.., AE, B”.. and BE, involving the source and 
vorticity distributions described in Section C of Chapter II are still useful. These are 
indeed time-independent coefficients since they are functions of geometrical parameters 
which are fixed in our rigid airfoil. Additional influence coefficients involving the shed 
vorticity panel and the wake core vortices must be defined. These coefficients need to 
be computed in each time step since their positions vary relative to the airfoil-fixed 
coordinate system. For that matter, as will be made clear later, those influence 
coefficients involving the shed vorticity panel must also be computed in every iteration 
within each time step for the same reasoning. 

a. More A's and B's Influence Coefficients 

Following the notations used previously in steady llow'. we define, with the 
use of the k- subscript to denote time-dependency, additional influence coefficients 
involving uniformly distributed singularities of source and vorticity. They are the A's 
and B's coefficients : 

® + •' normal velocity component induced at the i^ panel control 

point by unit strength vorticity distribution on the shed vorticity panel at time 
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• : tangential velocity component induced at the i^ panel control 
point by unit strength vorticity distribution on the shed vorticity panel at time 

^k* 

• (^'’^n + lj^k : X- velocity component induced at the shed vorticity panel mid 

point by unit strength source distribution on the panel at time tj^, 

• : y- velocity component induced at the shed vorticity panel mid 
point by unit strength source distribution on the panel at time tj^. 

• (B*'^n + lj^k : x-velocity component induced at the shed vorticity panel mid 
point by unit strength vorticity distribution on the panel at time tj^. 

• : y- velocity component induced at the shed vonicity panel mid 
point by unit strength vorticity distribution on the panel at time t,^. 

• • x-velocity component induced at the center of the h^ core 
vortex by unit strength source distribution on the panel at time tj^. 

• ■ y*'^siocity component induced at the center of the h^^ core 
vortex by unit strength source distribution on the panel at time t^^. 

: x-velocity component induced at the center of the h^ core 
vonex by unit strength vorticity distribution on the panel at time 

• (B''^^j)j^ : y-velocity com.ponent induced at the center of the h^ core 

vortex by unit strength vorticity distribution on the panel at time t.^. 

• (®'\n+l^k : X-velocity component induced at the center of the h^ core 

vortex by unit strength vorticity distribution on the shed vorticity panel at time 

'k- 

• + : y- velocity component induced at the center of the h^ core 

vortex by unit strength vorticity distribution on the shed vorticity panel at time 

b. New Os Influence Coefficients 

The presence of discrete core vortices in the wake requires the definition of 
new influence coefficients involving point singularity. They are the C's coefficients in 
our familiar notations : 

• (^*^im)k • velocity component induced at the i^ panel control 

point by unit strength m‘^ core vortex at time tj^, 

• tangential velocity component induced at the i^ panel control 
point by unit strength m’“^ core vortex at time . 

^^^'n-i-lm'k ■ x-velocity component induced at the shed vorticity panel mid 

point by unit strength m'^ core vortex at time tj^. 

• (O' 1 , _). : y- velocity component induced at the shed vorticity panel mid 

point by unit strength m core vortex at time tj^. 

• ■ x-velocity component induced at the center of the h^ core 
vortex by unit strength m^^ core vortex at time tj^. 



36 



* component induced at the center of the h'^ core 

vortex by unit strength core vortex at time 

2. Computation of Time-Dependent Influence Coefficients 

(B'^j and (B^ are computed exactly the same way as B"- and B'-j. 

are computed using Equations 2.14 and 2.15 with subscript n-r 1 replacing j. Similarly, 
and are calculated using Equation 2.12 with 9j set to zero and 

subsript i appropriately replaced. Also (A^^ 4 .^j)j^ and are calculated using 

Equation 2.13 with 0. set to zero and subsript i appropriately replaced. We do the 
same for ^ and j)j^ using Equations 2.14 and 2.15 respectively and so 

on for all the rest of A's and B's coefficients. The C's coefficients wUl be computed 
with different expressions from those of A's and B's because they are the velocities 
induced by unit strength core vortex. It can be shown easily, from the geometry of 
Figure 3.3, that their expressions take on the following forms, 




k 



- (9„),)i 



(eqn 3.11) 



(C'„)k = 



_ sin[9j - (9J^)| 
2" (fta)k 



(eqn 3.12) 



where : 



<Wk = V((xmi-xj2 + (ymi-yj2] 
xmj = */2(Xj + Xj + j) 

ymj = ‘/2(yj + yj + ^) 

x^ = X coordinate of m* core vortex at time t, 

m K 

Vrr, ~ y coordinate of m*-^ core vortex at time t. 

ill fC 



0j = arctan( ^ 






(^m^k " arctan( 



- V- 
xm. - y 

I 



m 



m 
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By the same token. and are computed by Equation 3.11 

while , ^)j, and (C^hm^k computed by Equation 3.12 if 6j is set equal to zero 

and the subscript i appropriately replaced. 

D. NUMERICAL SOLUTION SCHEME 
1. The Flow Tangency Conditions 

The flow tangency conditions of Equation 3.5 can be written using the 
influence coefficients as follows, 



n 



I [ (qj), 1 + Y, y B". + [ (V . n. ], 

j=l j=l 



k-1 



-i-LM. (B" ^,). +yf(CT )(T ,-r )] = 0 . i=l,2 n (eqn 3.13) 

^ * vv^k ^ i,n T- l^k ^ ^ ^ m-1 m' ^ ^ ^ ' 

m = 1 



where (^^.£301^1 evaluated by Equation 3.7 at the i^ panel control point 

This equation, though it seems complex, says nothing more than summing to 
zero all the velocity contributions due to individual singularity. Substituting (Yw)k 
Equation 3.2, collecting like terms and rearranging the equation into. 



n 



I ( A'-ij (Qj), ] = '/,[ (i!\) (B 

i = l 



n 

i,n + 1 




, i=l,2 n 

m = 1 



(eqn 3.14) 



2. The Iterative Solution Procedure 

Equation 3.14 is arranged in this form because we intend to solve 

j= 1.2 n) in terms of 7;^. Expanding Equation 3.14 for i= 1.2 n results in the 

following matrix equation. 



[ A ] { q }^ = Yk { B }^ + { C 



(eqn 3.15) 
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Figure 3.3 



Mid Point ofi^ Panel 

(xmj.ynij) 




Influence CoefTicients due to Point Singularities. 
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where [ A ] is an n x n matrix whose elements are known constants. { B and { C 
are n X 1 column vectors w’hose elements are known only if the shed vorticity panel at 
tj^ is established; that is, if and are known, then w'e can calculate ail the 
influence coefhcients on the right-hand-side of Equation 3.14. We therefore make use 
of this idea to formulate our iterative solution procedure as follows : 

(1) Project the wake core vortices downstream according to the time step and the 
local resultant velocities at their respective centers with respect to an inertial 
coordinate system. 

(2) Compute the coordinates of these core vortices relative to the airfoil-fixed 
coordinate system due to its time-dependent motion. 

(3) Start iteration cycle for current time step by initially assuming some guess 
values of Aj^ and 0j^. We can use, except for the first time-step, values 
obtained at previous time step. Compute then the influence coefficients needed 
in Equation 3.14 or 3.15. 

(4) Obtain (qj)j^ in terms of by solving Equation 3.15 as a linear system with 
two right-hand-sides by the same Gaussian elimination algorithm used in 
steady flow. 

(5) Calculate the tangential velocities at the trailing edge panels, all in terms of 

(6) Invoke the Kutta condition of Equation 3.6 (with some effons in algebriac 
manipulation) to solve for since it is the only unknowm in that equation. 

(7) Once is solved, (q^j^ are then known. We can then calculate the velocity 

components and (V^)j^ at the mid point of the shed vorticity panel. 

(8) Equation 3.3 and 3.4 hence enable us to update the values of Aj^ and 0j^. 

(9) Repeat the iteration cycle from steps (3) to (9) until converged values of Aj^ 

and are obtained. Alternatively convergence can be set for and 

(V^),. instead. 

(10) Compute the tangential velocities and disturbance potential at all panel 
control points in order to determine the pressure distribution which can be 
integrated to give forces and moments. 

(11) Compute the resultant velocities which occur at the centers of all the core 
vortices that will be convected down-stream. These velocities must be the 
absolute velocities with respect to an inertial coordinate system. 

3. Computation of Velocities 

The iterative procedure mentioned in the previous subsection requires 
calculation of tangential velocities at the trailing edge panels and the absolute velocity 
components (U^y)j^ and (V^^)j^. They are computed differently due to the use of a 
different frame of reference. 
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a. Tangential Velocities on Airfoil Panels 

The tangential velocities (i= 1,2,. ...,n) at ail the panel control points 

are calculated using the airfoil-fixed coordinate frame of reference as follows : 



+ , i=i,2 ,n (eqn3.16) 

m= 1 

b. Core Vortices Convection Velocities 

The resultani velocities at ail core vortices are calculated using the inertial 
frame of reference fixed with respect to but resolving them into components in the 
directions coincide with the airfoil-fixed coordinate system as shown below : 



a 



a 





n 



a 



(^h)ic = Z [ (A\j)k (qj\ 1 + Yk S 



j=i 




k-l 






(eqn 3.17) 



m = 1 
h 



n 



n 



(Vh), = E [ (A\p, (qp, 1 + Yk E 



j = l 



1 = 1 




k-l 






(eqn 3.18) 



m= 1 



h 



41 



Notice the use of instead iri Equations 3.17 and 3.18. Also 

the subscript h is just an index usable for any core vortex. We can obtain and 

if h is replaced by n-r I in these equations. 

4. Disturbance Potential and Pressure Distribution 
a. Why We Need the Disturbance Potential 

The concept of disturbance potential (p has been instrumental in the 
formulation of both the steady and unsteady flow problems. However, it has never 
gone beyond using it merely as a vehicle to understanding the superposition of simple 
flows. The disturbance potential need not be solved for at all in the steady flow 
problem formulation. This is because what one really is going after is the spatial 
derivative of this disturbance potential, i.e. the disturbance induced velocity, from 
which the pressure distribution can be obtained. We have, in all our solutions so far, 
been successful in avoiding any disturbance potential calculation since the concept of 
influence coefficients allows us a direct evaluation of the velocity. Unfortunately, as 
can be seen in Equation 3.10, when we proceed further to compute the pressure 
distribution on the airfoil surface in unsteady flow, we are faced with the oroblem of 
evaluating the disturbance potential (O. or more precisely the rate of change of (p, which 
we approximate by using a backward finite difference expression. Therefore, the 
pressure coefficients at the i"^ panel control point can be rewTitten, in terms of non- 
dimensional variables, as, 



KCJil, = - [(V‘).l 



2 _ «Pi)k-(<Pi)k-l 



h v-1 



(eqn 3.19) 



where (V^). is calculated by Equation 3.16 and is the non-dimensional (by 

Vqo) form of Equation 3.7 evaluated at the i^^ panel control point. 

We thus need to calculate at each time step, the disturbance potential at all 
the panel control points. Short of having to solve the Laplace's equation by a finite 
difference scheme, wc evaluate the disturbance potential ^p by integrating ‘he velocity 
field in two stages from upstream at infinity to the airfoil leading edge, then along the 
airfoil surface from the leading edge to each panel control point. Care must be taken 
here to include only the velocity contribution due to disturbances. 

One important question arises, in this approach, as to w’hat value of 
disturbance potential we should use at infinity before we carry out the line integral. We 
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must therefore analyse the behaviour of (p at iniiniiy by examining the singularities that 
constitute the disturbance. They are the source and vorticity distributions on the airfoil 
surface and the core vortices in the wake. These singularities induce no velocity at 
infinity from the knowledge of simple Qows. In other words, the disturbance potential 
(p at infinity is independent of spatial coordinates. The next question \ve should ask is 
whether (p at infinity is time-dependent? Let us adopt the view-point that if we are at 
infinity looking at our airfoil and its associated wake, we simply see a point vortex with 
a total circulation Fq at time tg. We have already identified that Fq remains constant 
by Helmholtz's theorem. It only gets redistributed, as time progresses, over the airfoil 
surface and in the wake. Notice that the previous statement regarding what one would 
see at infinity said nothing about the source distributions. The source distributions 
though vary (or get redistributed ) as the time progresses, the total source strength 
necessarily remains zero at all time in order to enforce a closed contour representing 
the airfoil thickness. This is also the reason why the unsteady flow solution needs an 
additional model to handle the vorticity conservation since the source conserv'ation is 
already implicitly so for a closed contour to exist. From these discussions, we are 
certain that the disturbance potential (p at infinity is an absolute constant (independent 
of time and spatial coordinates) whose value is fixed only by the initial condition we 
decide to start solving the unsteady problem. The actual value of (p at infinity is in fact 
immaterial so long as we know it is constant because its value disappears conveniently 
as we subtract (<Pi)j^.j from ((pj)j^ in Equation 3.19. 
b. Computation of Disturbance Potential 

We begin by choosing an arbitrary straight line extending upstream to 
infinity from the leading edge of the airfoil along a direction parallel to V^q. For 
practical purposes, we set infinity at say ten chord lengths away from the leading edge 
since the velocities induced, at field points thereafter, by the disturbances are small 
enough to be negligible. This line is divided into z panels with element lengths near the 
leading edge comparable to the panel sizes used on the airfoil. However, the panel size 
is progressively increased to take advantage of the inversely decaying induced velocities 
at larger distances. We compute the tangential components of the induced velocities at 
the mid points of these panels using influence coefficients analogous to those used on 
the airfoil panels. Using subscript f to denote these panel mid-points, we can define 
influence coefficients (A^^)j^, ^),^, and and compute 

them using the same expressions for calculating the A's, B's and C's coefficients used 



43 



before with cos0j replaced by ( — cosa), sin0; replaced by ( — sina) and subscript i 
replaced by f. With the help of these influence coefficients, the tangential velocity 
induced by disturbances at the panel mid point is : 



n 



= I ( 1 " •/]< 1 

i=l i=l 



k-1 



+ (y.A H- 1 ( (cvj, (r^.j - r^) i 

m = I 



(eqn 3.20) 



valid for f= The disturbance potential at the airfoil leading edge is the sum of 

the products of the disturbance induced velocity at each panel and the panel length. 



(<f>iA = - S [(''Vr'k f 1 “ ^ 1 “ M-)* I 



(eqn 3.21) 



f= 1 



Similarly, for the line integral over the airfoil surface, we compute the 
tangential component of the disturbance induced velocity at the i^ panel control point 
using the following equation : 

t(V<p);)k = S t A'ij (qj), 1 + T, i B>.j 

i=l j = l 

+ (V J, + ,), + !:[ (C‘ (r„., - r„) 1 (eqn 3.22) 

m — 1 

which is valid for i= l,2,....,n. Performing the line integration by summation, the 
disturbance potential at the i^ nodal point on the airfoil is : 

for a > i ^ ij^ 

for i|g > i > 1 (eqn 3.23) 



^^nodc i\ ~ ~ "j,j + 1 

j = 'le 
ile-' 



(<p,A ' 



J 



where rj j + ^ denotes the panel length. 
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rjj + , = Vl(Xj + ,-Xj)2 + (Vj + j-yp^i 



Finally, the disturbance potential at the i*** panel control point is, 



(<Pi)„ = [ (<P„<,de i>k + ('Pnod=i+l)kl 



1 1,2, ....,11 



(eqn 3.24) 



5. Computation of Forces and Moments 

The Cg, and about the leading edge are calculated in exactly the same 
way as it is done for the steady flow problem by integrating the pressure distribution 
(See section D-4 of Chapter II). 

E. FLOW MODELING OF SHARP EDGE GUST FIELD 

The unsteady flow solution described so far can be extended to the study of 
airfoils penetrating a sharp edge gust by modifying the boundary conditions with the 
assumption that the gust front remains straight while passing through the airfoil. The 
same assumption has been used in both [Ref. 3] and [Ref. 4]. An additional model in 
[Ref 3] using distribution of singularities along the gust front had successfully 
attempted to simulate the distortion of the gust front passing over the airfoil surface. 
It was shown that the pressure distributions, during the time when the gust front 
remained on the airfoil surface, were affected only at the neighbourhood of the gust 
front. The overall pressure upstream and do'v^mstream of the gust front stayed 
essentially the same. The distorted gust front model is not used in program U2DIIF. 
The use of the relatively simple yet sufficiently accurate model of a straight gust front 
affords the modifications to the unsteady flow solution to be confined only to the flow 
tangency conditions. That is to say, the expression of in Equation 3.7 would 

include the gust velocity for panels that are already in the gust field during the 
penetration phase. Similarly, the computation of core vortex velocities using 
Equations 3.17 and 3. IS have the gust velocity added to the if the core vortices are 
already in the gust held. 

In an attempt to generalise the solution for cases where airfoils enter the gust 
field at an angle of attack, the convenient model used in [Ref 3] by setting the 
computation to proceed, for the undistorted gust front simulation, so that the gust 
front always coincides with the nodal points is difficult to implement. At any one time. 
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if an airfoil enters a gust field at an angle of attack, the gust front would appear in 
between two nodes of a particular panel on one surface while the gust front proceeds 
from node to node on the other surface. We therefore further modify the flow 
tangency condition only on that particular panel where the gust front lies in between 
two nodes by taking the gust velocity on that panel to be proportional to the fraction 
of panel length partially submerged in the gust field. 
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IV. DESCRIPTION OF COMPUTER CODE U2DIIF 



A. PROGRAM U2DIIF STRUCTURES AiSD CAPABILITIES 
1. Restrictions and Limitations 

% 

The numerical formulations of both the steady and unsteady flow problems 
outlined in the previous Chapters are coded in a FORTRAN computer program called 
U2DIIF (See Appendix A for the program listings). The present solution methods 
treat the inviscid and incompressible flow as an approximation to the real flow so long 
as the viscous effect is negligible and the flow stays attached on the airfoil surface at all 
time. These restrictions are no strangers to any one who is familiar with any other 
potential flow solution methods. Other than the implicit restrictions of potential flow 
solution, the method is entirely general in that the shape of airfoil is arbitrary and any 
arbitrary continuous motion of the airfoil could be simulated using either the closed 
form (i.e. explicit equations) or discrete data points to describe the time-histor}/’ of the 
translational and rotational velocities. 

The storage of the computer that carries out the calculations may be the other 
limitation one should consider. The storage requirements grow rapidly with the number 
of panels (n) and the number of computation time steps (m). By far the prime 
contributor to this storage requirement comes from the massive amount of influence 
coefficients. The number of influence coefficients increases with the square of the 
number of panels (n^). Each time step increment adds (2n + m^) more influence 
coefficients due to the formation of shed vorticity. The current program fixes the 
maximum number of airfoil panels to 200 and the maximum allowable time steps is 
also 200. 

An additional constraint worth mentioning concerns the gust field simulation 
whereby the current solution methodology is valid except in the use of the same 
pressure equation arising from Che unsteady Bernoulli's equation (See Equation 3.10). 
The fundamental assumption underlying the derivation of this equation is the 
irrotationality of the flow field. There is no doubt that the flow fields upstream and 
downstream of the gust front are irrotational. However, when one needs to obtain the 
pressure on the airfoil surface, an implicit integration is done across the gust front. A 
flow field inclusive of the gust front is rotational since the line integral of velocity in a 
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closed path does not vanish when the gust front is crossed. Failing the proper 
derivation of a new pressure equation applicable to unsteady rotational flows, care 
must be exercised to regard the present method as an approximate solution to gust 
fields of weak strengths only. 

2. Current Structures of U2DIIF MAIN Program 

The overall program logic-flow chart is as shown in Figure 4.1. The program 
first reads in the input data from filecode 1 and sets up the airfoil panel nodes and 
slopes. Immediately after that, the steady flow calculations are executed for the initial 
angle of attack Uj according to the solution scheme described in Section D of 
Chapter II. The steady flow solution is included primanly to : 

• Provide the necessary initial parameters for the unsteady flow solution to 
proceed in time. In other words, the steady flow solution handles the and 
initial angle of attack one decides to begin the unsteady flow calculation. 

• Allow the code to function directly as a steady How solver as and \vhen 
necessary without having to do the time consuming unsteady flow iterative 
solution and approach the steady flow as time approaches infinity. 

The program terminates once the. steady flow calculations are done if the 
program determines, based on the input data set by user, no requirement for unsteady 
flow solutions. Otherwise the unsteady flow calculations will be activated by selecting 
and computing the rigid body motions of the airfoil and the corresponding 
computation time-step size. Currently, all the time dependent motions are equation- 
generated, they are the positions and rates of the translational and rotational motions. 
Incorporated as case-runs within the program U2DIIF are the following motions : 

(1) Step change in angle of attack from any initial value. 

(2) Modified-ramp change in angle of attack about any pivot point from any 
initial value. 

(3) Harmonic translational motion at any angle of attack. 

(4) Harmonic rotational motion about any pivot point at any mean angle of 
attack. 

(5) Sharp edge gust penetration at any angle of attack. 

Should one decide to generate the airfoil's motion using discrete data points as a 
function of time, the program could be easily modified. 

The computation time-step sizes for the harmonic translational and rotational 
motions are constant values determined by the frequencies (FREQ) and the number of 
computation per cycle (DTS). For the case of step change in angle of attack, the 
computation time-step size is progressively increasing, from a starting value (DTS), as 
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Figure 4.1 Flow Chart for U2DIIF Computer Code. 
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Figure 4.1 .(corn'd.) 
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time increases. The case of the modified ramp change in angle of attack, adopts 
initially a constant computation time-step size (DTS) during the transient rising of the 
angle of attack. Once the final angle of attack is reached, the computation time-step 
size is progressively increased also. Similarly for the case of airfoil penetrating a gust 
field, the computation time-step size (DTS) is constant during the period when the gust 
front remains on the airfoil surface but progressively increases once the entire airfoil is 
submerged in the gust field. These variations in computation time steps are to provide 
greater flexibility both in capturing transients and covering relatively large total time of 
computation without having to contend with the storage space requirements described 
previously. These variations in time-step sizes described so far are associated with 
setting the input parameter TADJ to zero. If TADJ is chosen to be non-zero, all the 
case-runs would compute initially using the starting time-step sizes, based on DTS for 
non-oscillatory motions and FREQ & DTS for harmonic motions, and the program 
would prompt for an user choice of time step adjustment. If the answer is yes, the 
program would back-track the previous solution and recompute the current solution 
using an adjusted time-step size that is TADJ times the initial value (DTSX- This special 
time step variation feature gives the program added capability of allowing an 
interactive time step selection during the progress of unsteady flow computation. The 
ability to back-track and recompute the current solution using a different time-step size 
enhances the possibility of using program U2DIIF together with a viscous flow' solver 
forming an Inviscid- Viscous-Interactive solution scheme which often requires such time 
step variations. 

The MAIN program performs the iterative solution procedures set out in 
Section D of Chapter III. The convergence check during the iterative solution is done 
through the user specified tolerance between successive iterative solutions of both 
(Uw)j^ and (V^)j^. The solution continues into the next time-step by selecting the time 
step size according to the particular case-run and projecting all the wake core vortices 
downstream so that their new positions relative to the airfoil at the new' time step can 
be correctly determined. 

B. DESCRIPTION OF SUBROUTINES 
1. Subroutine BODY 

This subroutine is called by subroutine SETUP if the user selects an airfoil 
that is either a NACA XXXX or 230XX type. It in turn calls subroutine NACA45 to 
obtain the airfoil thickness and camber distributions and returns with the computed 
(x,y) coordinates of the panel nodal points. 
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2. Subroutine COEF 

This subroutine is called by the MAIN program in the unsteady flow 
calculations. It utilises, at each iteration cycle, the influence coefficients generated by 
subroutine INFL to calculate the coefficients of the matrix Equation 3.15 by expanding 
Equation 3.14. These matrix coefficients are necessarily set up in this way so chat the 
source strengths could be solved in terms of the vorticity strength by subroutine 
GAUSS as a linear system with two right-hand-sides. 

3. Subroutine COFISH 

This subroutine is called by the MAIN program to set up the coefficients of 
the matrix system of Equation 2.18 for steady flow where the source strengths and 
vorticity strength are solved simultaneously by subroutine GAUSS as a linear system 
with one right-hand-side. The matrix coefficients are calculated using Equations 2.16 
and 2.17. 

4. Subroutine CORVOR 

This subroutine is called by the MAIN program at nearing the end of the 
unsteadv flow calculations before stanins a new time steo. It comoutes the resultant 
convective velocities for all the wake core vortices with respect to an inertial frame of 
reference according to Equations 3.17 and 3.18 where all the appropriate influence 
coefficients are linked through common block from subroutine INFL. 

5. Subroutine FANDM 

This subroutine is used in both the steady and unsteady flow calculations. It 
is called by the MAIN program immediately after the pressure distribution over the 
airfoil panels are known so that it can perform the simple integration of pressure in the 
appropriate directions to give the aerodynamic force and moment coefficients of lift, 
drag and pitching moment about the leading edge according to Equations 2.22 through 
2.26. 

6. Subroutine GAUSS 

This subroutine is the standard linear system solver that employs the well 
known Gaussian elimination with partial pivoting and operates simultaneously on a 
user specified number of right-hand- sides. It is called by the MAIN program m both 
the steady and unsteady flow calculations. In order to use GAUSS, the coefficients of 
the augmented matrix must be set up so that GAUSS will return the solutions 
replacing the corresponding columns of the augmented matrix that were initially 
occupied by the right-hand-sides. The coefficient set-ups are done by subroutines 
COFISH and COEFF respectively for the steady and unsteady flow problems. 
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7. Subroutine INDATA 

This subroutine is called by the MAIN program to read in the first three sets 
of data cards and returns to the MAIN program if IFLAG * 0. Otherwise it continues 
to read in the fourth data card as the NACA number corresponding to the type of 
airfoil and calculates the thickness parameters that will be used by subroutine 
NACA45. 

8. Subroutine INFL 

This subroutine is the generator for all the influence coefficients that need to 
be stored and used by many subroutines associated with the unsteady flow calculations. 
It utilises the known relative geometrical parameters of the singularities to carry out 
computations based on Equations 2.12 through 2.15, 3.11 and 3.12. The MAIN 
program calls this subroutine in every iteration cycle of each time step so that the time- 
dependent coefficients can be updated as and when necessary. Time-independent 
coefficients are computed only once in the entire unsteady flow solutions. Those 
influence coefficients involving the wake core vortices are updated in each time step 
while those involving the shed vorticity panel are calculated as frequently as the 
number of iterations take to terminate a converged solution. It, however, does not 
compute and store those influence coefficients needed for the determination of 
disturbance potential (Equation 3.20) simply because they are used only once in each 
time step. 

9. Subroutine KUTTA 

This subroutine is called, in the unsteady flow calculations, by the MAIN 
program during every iteration cycle in each time step to invoke the Kutta condition 
for unsteady flow expressed in Equation 3.6. It calculates the tangential velocities at 
the trailing edge panels using Equation 3.16 in terms of the unknown vorticity strength 
that is manipulated and solved algebraicly. 

10. Subroutine NACA45 

This subroutine is called bv subroutine BODY if the airfoil selected belonas to 

✓ w 

the families of NACA 4-digits airfoils or the NACA 5-digits airfoils of type 230XX 
who share common thickness distributions with the 4-digits airfoils having the same 
thickness to chord ratio. The thickness and camber distribution data of these airfoils 
are calculated and returned to BODY. 
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11. Subroutine PRESS 

This subroutine is called by the MAIN program to calculate the pressure 
distribution over the airfoil panels after the iterative solution for the unsteady flow 
problem has successfully met the convergence criterion. It first computes the 
tangential velocities at all panel control points using Equation 3.16, then performs the 
disturbance potential evaluation at the current time step according to Equations 3.20 
through 3.24. Together with the disturbance potential data obtained from the previous 
time step, it calculates the pressure distribution using Equation 3.19. 

12. Subroutine SETUP 

This subroutine sets up the panel nodal coordinates for MAIN program by 
reading the 4^*^ and 5^ data sets of the input file if IFLAG = 1 is set. It skips the data 
reading if IFLAG = 0 and proceeds to set up the node distribution and call subroutine 
BODY to calculate the airfoil coordinates. The node distribution adopts a cosine 
formula in order to have closelv oacked oanels towards the leadins and trailina edges 
for improvements in solution accuracy. Regardless of how the nodal coordinates are 
obtained, SETUP determines the panel slopes and airfoil perimeter length. 

13. Subroutine TEWAK 

This subroutine is called by the MAIN program at every iteration cycle of 
each time step of the unsteady flow calculations to compute the resultant velocity 
components at the mid point of the shed vorticity panel using Equation 3.17 and 3.18 
These velocity components are necessary to ensure the correct establishment of the 
shed vorticity panel length and orientation which governs the successful 
implementation of the iterative solution scheme for the unsteady flow problems. 

14. Subroutine VELDIS 

This subroutine returns to the MAIN program the velocities and pressure 
distributions for steady flow calculation using Equations 2.20 and 2.21. It also 
performs the evaluation of the disturbance potential at the panel control points. 
Though these disturbance potential data are not necessary for steady flow solution, 
they will be needed in the first time seep of the unsteady How pressure calculation. 

C. INPUT DATA FOR PROGRAM U2DIIF 

Program U2DIIF reads its input data from filecodc 1. An example of the input 
data file is as shown in Appendix B for the case where the airfoil nodal coordinates arc 
input by user. User could however let the program generate the nodal coordinates if 
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the airfoil chosen happens to belong to the family of NACA 4-digits or 5-digits of type 
230XX. To do this, simply change I FLAG to zero in the first item of the 3'’^^ set of 
data card and replace the nodal coordinates data in the 4^ and sets of data cards 
by a single data card containing only the particular airfoil's NACA number using 
Format (15). Figure 4.2 contains an itemised description of the sequential input 
variables. 

D. OUTPUT DATA FROM PROGRAM U2DIIF 

Appendix C contains a sample output data generated by using the input data set 
shown in Appendix B. Due to the repetitive nature of output as the computation time 
progresses, only data at selective time steps are shown. The output data file begins with 
writing out what the program has read from the input data file followed by the 
computed nodal coordinates only if they are program generated, otherwise proceeds to 
write the computed airfoil perimeter length. The next set of output data are the steady 
flow solution parameters of distributed source strengths, vorticity strength, pressure 
and velocity distributions as well as the force and moment coefficients. The output data 
terminates at this point unless unsteady flow solution is required. It then prints, for 
each time step, the unsteady flow solution parameters similar to the previous output 
for steady flow with additional information pertaining to the rigid body motions and 
trailing wake vortices data. An explanation of the output variable names are listed in 
Figure 4.3. All output parameters are non-dimensional quantities. 
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Data Set #l 
ITITLE 

Data Set f^2 
TITLE 

Data Set #3 
I FLAG 

NLOWER 

NUPPER 

Data Set #4 
X(I) 

Data Set §5 
Y(I) 

Data Set #6 
ALPI 
DALP 



TCON 

FREQ 

PIVOT 

UGUST 

VGUST 

Data Set #7 
DELHX 

DELHY 

PHASE 

Data Set US 
TF 
DTS 



TOL 

TADJ 



Format (15) - I data card 

- Number of title cards to be used in Data Set =2. 

Format (20A4) - ITITLE data cards 

- Headings to be printed on output for case run identification. 

Format (315) - 1 data card 

- 0 if airfoil is NACA XXXX or 230XX. 

- I otherwise. 

- Number of panels used on airfoil lower surface. 

- Number of panels used on airfoil upper surface (need not 
be the same as NLOWER). 

Format (6F10.6) if IFLAG= 1 - variable data cards 

- x-nodal coordinates (divided by the chord length, c). A total 
of n + I nodal points divided irito 6 points per data card. 

Format (6F10.6) - variable data cards. 

- v-nodal coordinates (divided bv c) correspondina to the 
Data Set #4 if I FLAG = 1. 

Format(7F10.6) - 1 data card 

- Initial angle of attack (AOA) in deg. 

- Increment in AOA in dec for non-oscillaton.* motions. 

- Maximum amplitude of AO.A. change m deg' for rotational 
harmonic motions. 

- Non-dimensional rise time (V^qI/c) of AOA for motion 
involving modified-ramp change in AOA. 

- Non-dimensional oscillation frequency (o)c/Vqq) for 
harmonic motions. 

- The length from the pivot point to the leading edge divided 
by c (postive aft) for rotational motions. 

- Magnitude of gust velocity (divided by Vqq) along Vqq- 

- Magnitude of gust velocity (divided by Vqq) perpendicular to Vqq 

Format (3F10.3) - 1 data card 

- Amplitude of chordwise translational oscillation divided by c. 
(positive forward). 

- Amplitude of transverse translational oscillation divided by c. 
(positive downward). 

- Phase angle in des between the chordwise and transverse 
translational oscillation with the latter as reference. 

Format (4F10.3) - 1 data card 

- Final non-dimensional time to terminate unsteadv (low solution. 

- Startins time-steo size for non-osciilatorv motions if T.ADJ = 0.0. 

- .Numbe'r of computation steos per cycle for harmonic motions. 

- Baseline time-step size for all motions if TADJ 5=0.0. 

- Tolerance criterion for checking the convergence between 
successive iterations of (U ). and (V ), 

W 1C 1C 

- Factor by which DTS will be adjusted. 



Figure 4.2 List of Input Variables. 



56 



TK 


- Time step tj^. 


TKMl 


- Time step 


ALPHA(T) 


- Angle of attack at time tj^. 


OMEGA(T) 


- Rotational velocity (positive counter clockwise) at time tj^. 


U(T) 


- Chordwise translation velocity (postive forward) at time tj^. 


V(T) 


- Transverse translational velocity (positive downward) at time t^. 


NITR 


- Iteration number. 


vocw 


- Iterative solution of (U^^.)j^. 


VYW 


- Iterative solution of 


Wake 


- Iterative solution of shed vorticity panel length Aj^. 


THETA. 


- Iterative solution of shed vorticity panel orientation 0,^. 


GAM. VI A 


- Iterative solution of the strength of vorticity distribution. 


J 


- Panel number. 


X(J) 


- x-coordinate of the mid point of panel. 


Q(J) 


- Strength of source distribution on the panel. 


CP(J) 


- Pressure coefficient at the mid point of panel. 


V(J) 


- Total tangential velocity at the mid point of j* panel 
referenced to the airfoil-fixed coordinate system. 


CD 


- Drag coefficient. 


CL 


- Lift coefficient. 


CM 


- Pitching moment coefficient about the leading edge. 


M 


- Trailing wake core vortex number. 


X(M) 


- x-coordinate of the center of m^ core vortex. 


Y(M) 


- y-coordinate of the center of m^ core vortex. 


CIRC 


- Circulation strength of the m* core vortex. 



Figure 4.3 List of Output Variables. 
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V. RESULTS AND DISCUSSIONS ON CASE-RUNS 



This Chapter presents the results of numerous case-runs of U2DIIF code for the 
purpose of verifying the code. The various airfoils used in the case-runs are deliberately 
chosen to be the same as those airfoils where direct comparison ofresults can be made 
with either theoretical analyses and/or experimental data available in the literature. 

A. STEP CHANGE IN .ANGLE-OF-ATTACK 

1. Case-Run Definitions 

Consider an airfoil initially at zero angle of attack to the free stream Voo that 
undergoes a step change in angle of attack at tg. The resulting flow should then 
provide the time-dependent information on the build-up of aerodynamic forces and 
moments on the airfoil resembling the classical results of Wagner [Ref. 5] calculated 
based on linearised theory. Although Wagner prescribed a slightly different initial 
condition in that the airfou is initially at rest and impulsively started at an angle of 
attack and velocity Vqq, the difference is insignificant, especially for a symmetrical 
airfoil. This is because the seemingly different initial conditions when translated into 
the mathematical model means that the step change in AOA uses non zero initial 
circulation Tq at tg with non zero initial disturbance potential at infinity if the airfoil is 
cambered. For a symmetric airfoil, these initial values are all zeroes and therefore 
mathematically would be the same as the initial conditions prescribed by Wagner. 

2. Results and Discussions 

a. Von Mises 8.4% Thick Symmetrical Airfoil 

A 8.4% thick symmetrical Von Mises airfoil is used for this case-run where 
the airfoil performs a 0.1 rad (or 5.73 *’) step change in AOA. Figure 5.1 illustrates the 
changes in the pressure distributions over the airfoil at time instances corresponding to 
the airfoil having traveled distances, in terms of chord length, of 0.2. 0.5, l.O, 2.0 and 
The associated trailing wake patterns at these time instances Tess t= are shown 
in Figure 5.2. The time-dependent build-up of aerodynamic coefficients of lift, drag, 
pitching moment and the circulation strength over a computation period of two 
traveled chord length are shown in Figure 5.3. Notice that the lift, pitching moment 
and circulation results are normalised by the respective steady state values at the same 
AOA. The apparently large initial loading on the airfoil shown in Figure 5.3 correlates 
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consistently with the results of the Piston Theory of [Ref. 6] which predicts the starting 
load on an arbitrary wing to be 

Ct = 4a/M 

^starting 

where M is the Mach number. In the case of an incompressible flow (M = 0), the initial 
loading would be infinitely large. The same large initial loading was obtained by Kim 
and Mook [Ref. 7] who used continuous vorticities as panel singularities instead of our 
source and vorticity approach. Perhaps what remains most surprising is that the work 
of Basu and Hancock [Ref. 3] did not predict this initial loading, although they used 
the same singularity distributions as U2DIIF code. The initial large loading in lift and 
pitching moment decreases rapidly over a short time span, whereby the airfoil traveled 
approximately one-tenth chord length, before rising in a manner parallel to the Wagner 
Function. The drag, however, continues to decrease monotonically after the initial 
sharp fall. The circulation rises, as continuous shedding of vorticity takes place, slowly 
from the initial condition of zero to the asymptotic steady state value as time 
approaches co. These results, disregarding the initial large loading associated with 
incompressible flow, are in close agreement with the results of [Refs. 3,4,7]. 
b. Thickness Effects on the Wagner Function 

In order to more closely correlate the results of U2DIIF code to the 
theoretical prediction of Wagner, we performed the step AOA change calculations for a 
very thin (1% thickness) NACA 4-digit symmetrical airfoil which in reality should 
represent a flat plate. The results are plotted as shown in Figure 5.4. Shown also on the 
same Figure are the results of the 8.4% thick Von Mises airfoil and a 25.5% thick 
symmetric Joukowski airfoil. The initial loading falls off less rapidly for the case of the 
simulated flat plate as compared to other thick airfoils but the subsequent rise in lift 
follows very closely the Wagner Function. 
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Figure 5.1 Pressure Distributions at V'arious Time Instances 
Resulting from a 0.1 rad Step Change in AO.A for a 
S.4'^/0 Thick Synunctric \'on Mises .Airtbil. 
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(b) tVoo/c = 0.5 



Figure 5.1 . (cont'd.) 
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Figure 5.1 . (cont'd.) 
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(d) tVoo/c = 2.0 



Figure 5.1 . (cont'd.) 
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Figure 5.1 . (conc'd.) 
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Figure 5.2 Trailing Wake Patterns at Various Time Instances 
Resulting from a 0. 1 rad Step Change in AO.A for a 
8.4% Thick Symmetric Von .Mises Airfoil. 
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Fi'jurc 5.2 . (corn'd.) 
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(b) Normalised Pitching Moment vs iV^;c. 
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Figure 5.3 
Resulting 
8.4% 



Time-Dependent Aerodynamic Parameters 
from a 0.1 rad Step Change in AOA for a 
Thick Symmetric Von Mises Airfoil. 
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(d) Normalised Circulation r/Too 



Figure 5.3 . (cont'd.) 
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Vy'’agner Function 
1% Thick NACA-0001 Airfoil 
S.4% Thick Von Mises Airfoil 
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Figure 5.4 Time-Dependent Lift Resulting from 
Step Change in AOA for Airfoils 
of Various Thicknesses. 
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B. MODIFIED RAMP CHANGE IN ANGLE-OF- ATTACK 



1. Case- Run Definitions 

The case of a step change in AOA can be considered as an useful check for 
U2DIIF code since a handful of results from other theoretical analyses are available. 
However, a step change in AOA is practically not realisable since all motions, short of 
having infinite velocities, take place over a finite time span. The step change in AOA 
that is practically possible is in fact some form of ramp rise over a short time span with 
large velocity. Even so, due to the inertia of the airfoil, an exact ramp rise in AOA 
does not physically describe the actual motion of the airfoil since finite time is also 
involved before the airfoil could build up its ramp velocity. Same argument holds at the 
end of the ramp rise before the airfoil could stop at the final value of AOA. Therefore 
a so called modified ramp, with some form of rounding at the two ends of a ramp, is 
more likely to describe anything close to what is physically achievable. The theoretical 
work of Homentcovschi in [Ref 8] considers the case of a flat plate that moves with 
constant velocity and changes the incidence about the mid chord, through a particular 
ramp fashion, described mathematically as. 



a(t) = 



'0 for t < 0, 

J 5a (3 - 2t/T) t-/x^ 
j5a for t > T 



for 0 ^ t ^ I 



where 6a is the magnitude of the AOA change and T is the rise time for the AOA to 
reach its final value. This particular function, plotted as shown in Figure 5.5, does in 
fact describe such a modified ramp. 

2. Results and Discussions 
a. Flat-Plate Case-Run 

Since the results of [Ref 8] serves as another excellent source for the 
verification of U2DIIF code, the obvious thing to do is to use U2DIIF to compute for 
the case of a fiat plate, again simulated by the l°'o thick 'N'AC.\-000l airfoil, e.xecuting 
this modified ramp rise of 0. 1 rad .AOA over a rise time of l.i chord length. This rise 
time is chosen simply to facilitate a direct comparison of results to [Ref 8] which used 
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a rise time of 3 half-chord lengths. The results of computation are shown in Figure 5.6 
and 5.7. Figure 5.6(a) takes a close look at the build up of lift during the initial period 
when the airfoil moves a distance of six chord lengths. The lift initially rises to about 
82% and then decreases to about 66% of the steady state value during the transient 
rise time. Thereafter it increases monotonically in a manner parallel to the Wagner 
Function. Figure 5.6(b) is a zoom view of the rather slow convergence of lift to the 
steady state value. It takes the airfoU to cover a distance of around 50 chord length 
before the lift builds up to almost 99% of the steady state value. The same results were 
obtained in the theoretical analysis of [Ref 8]. Figure 5.7 shows a collection of the 
time-dependent aerodynamic parameters resulting from this particular case-run. 
b. Thickness Effects 

The same modified ramp function is used on the 8.4% thick Von Mises 
airfoil. The resulting lift-history plotted in Figure 5.8 shows a lower peak value of lift 
during the transient AOA rise as compared to the case of a Hate plate though a similar 
trend of lift rise is obtained. Figures 5.9 and 5.10 show the results of pressure 
distributions and trailing wake patterns at various time instances. One could directly 
compare these Figures to the corresponding Figures arising from the step AOA change 
calculations and see the remarkable differences in transient characteristics as a result of 
varying the prescribed motions. Incidentally, one should realise that the non- 
dimensional rise time of 1.5 chord length is a deceivingly large number. In fact, when 
one converts this to the real time for an airfoil of 10 ft chord length moving at a low 
Mach number of 0.2, the rise time is indeed only of the order of 0.06 sec. which for 
practical purpose is close enough to a step. Neverthless, the transient part of the lift 
response is entirely governed by how one prescribes the transient motion. 
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ANGLE OF ATTACK VS TIME 




Figure 5.5 The VI edified Ramp AOA Change. 





Figure 5.6 Normalised Lift Resulting from a 

Modified Ramp AOA (6a = 0.1 ra&, t = 1.5) about 
the .Mid Chord of a NACA-0001 Airfoil. 
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(a) Normalised Lift C^. C^ vs tV-o/c. 
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Figure 5.7 Time-Dependent Aerodynamic Parameters Resulting from a 
Modified Ramp ;\OA (5a = 0.1 rad. r = 1.5) about 
the Mid Chord of a NACA-0001 airfoil. 
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(c) Drag vs rV^Q/c. 




Figure 5.7 . (cont'd.) 
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Figure 5.8 Normalised Lift Cj. Cf Resulting from a 
.Modified Ramp AOA (6a = 0.1 r3^. r = 1.5) about 
the .Mid Cliord of a .Vlises S.4'?o Tliick airfoil. 
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(a) cV'oo/c = 0.2 



Figure 5.9 Pressure Distributions at Various Time Instances Resulting from a 
Modified Ramp AOA (6a = 0.1 rad, t = 1.5) about 
the Mid Chord of a Mises 8.4% Thick airfoil. 
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(b) tVoo/c = 0.5 



Figure 5.9 . (cont'd.) 
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Figure 5.9 . (cont'd.) 
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(d) tV^/c = 2.0 



Figure 5.9 . (cont'd.) 
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(b) cV3o,c = 0.5 VC 



Figure 5.10 Trailing Wake Patterns at Various Time Instances Resulting from a 
Modified Ramp AO A (6a = 0.1 rad, t = 1.5) about 
the Mid Chord of a Mises 8.4% Thick airfoil. 
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Figure 5.10 . (conc'd.) 
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c translational harmonic oscillation 



1. Case- Run Definitions 

Although the U2DIIF code is capable of computing unsteady flow solution 
for any general translational harmonic motion described by a chordwise and a 
transverse components bearing a given phase relationship, 

hy(t) = 6hy sin(cot) 
h^(t) = 6h^ sin((Ot + X) 

where (0 is the oscillation frequency, X is the phase angle between the two oscillation 

components and 6h & 6h are the magnitudes of chordwise and transverse oscillations 

X y 

respectively. The case-run to be presented in this section selects the motion to consist 
of only the transverse component, i.e. the heaving or plunging motion. A NACA-0015 
airfoil is chosen for the case-run. The airfoil is initially at zero AOA with the 
freestream Vqq and performs the plunging oscillation at an amplitude of 6hy and a 
reduced frequency of coc/Vqq 

2. Results and Discussions 

Figures 5.11 and 5.12 present the results of an airfoil executing a plunging 
motion at an amplitude of 0.018c but with two different reduced frequencies of 4.3 and 
17.0 respectively. These numbers are chosen to coincide with those numbers used in 
[Ref. 4]. Excellent correlations are obtained. Notice from these Figures that the 
oscillation frequency has a great influence on the magnitudes of the aerodynamic 
parameters due to the formation of significantly different trailing wake patterns for the 
same oscillation amplitude. Also to note is that the width of the resulting trailing w'ake 
is much larger than the amplitude of the oscillation, reinforcing the fact that the 
unsteady flow is strongly governed by the shed vorticity in the trailing wake. The lift 
and pitching moment oscillate at the same frequency as the airfoil motion but slightly 
out of phase, the phase differences vary with the oscillation frequency. The drag is 
however oscillating at about twice the frequency of the airfoil motion with a negative 
mean value, indicating that the plunging action indeed generates some propulsive 
thrust. The same conclusion was arrived at in the experimental work of Halfman 
[Ref 9] using a symmetrical NACA-0012 airfoil. 
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(a) Position of Airfoil (Positive Downward) 




I'oat) (Zz) 

(b) Time History of Cj, and over 2 Cycles. 



Figure 5.11 Harmonic Plunging .Vlotion of a NACA-00I5 Airfoil 
6h^. = 0.01 Sc. (oc.'V ^ = 4.3. 
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Figure 5.11 . (cont'd.) 
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(a) Position of Airfoil (Positive Downward) 




vO)t) ( Z/T) 

(b) Time History of Cf, and over 2 Cycles. 



Figure 5.12 Harmonic Plunging Motion of a N.ACA-0015 .Airfoil 
oh = O.OlSc, 0)c = 17.0. 
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Figure 5.12 . (cont'd.) 
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D. rotational harmonic oscillation 



1. Case- Run Definitions 

The case of an airfoil oscillating harmonically in pitch will be uniquely defined 
only if a pivot point is prescribed. If a fixed pivot point is used in the calculations, take 
for example the leading edge, any pitching motion about a pivot point other than the 
leading edge would need to be described as composed of a pitching and a translation 
motions about the leading edge. Program U2DIIF handles such conversion 
automatically without having the user to figure out the combined motion. This applies 
also to the case of modified ramp rise in AOA. The harmonic pitching oscillation is 
described by, 

a(t) = 5a sin(o)t) 

where 5a and (O are the amplitude and frequency of oscillation respectively. 

2. Results and Discussions 

The results for the case of the S*.4% thick Von Mises symmetric airfoil, 
oscillating at an amplitude of O.Ol rad (or 0.573°) at a rather high reduced frequency 
of (oc.'Vqq = 20.0 about the leading edge, are shown in Figure 5.13. The lift, drag and 
pitching moment time-history as well as the trailing "wake patterns are very’ much 
similar to the case of a plunging airfoil at frequency of the same order of magnitude. 
The differences are clearly in the magnitudes and phase angles. These results check 
closely to those of [Ref 3]. Figure 5.14 shows the results of the same Mises airfoil 
performing another harmonic oscillation at a lower reduced frequency of 0.8 and a 
large amplitude of 0.3973 rad about a pivot point 0.5 chord length ahead of the leading 
edge. [Ref 4] conducted the same analysis for this pitch oscillation although the 
reason for using such a high amplitude of almost 23° was not clear. It is envisaged 
that such high amplitude may result in fiow separation. Nevertheless, the case-run is 
carried out assuming validity of attached flow for the sake of comparing the results. 
Perhaps an inherent advantage, in the light of U2D11F code verification, with the use 
of high amplitude in this case-run is that a discrepancy, if any, would show up 
significantly. Due to the use of difierent computation time-step size the pressure 
distributions on the airfoil, shown in Figure 5.14, do not correspond one-to-one at 
exactly the same angular positions as those presented in [Ref 4]. However, the angular 
positions are matched to within 0.001°, 0.05° and 0.8° respectively for the three 
pressure distributions shown. They correlate very well to the results of [Ref 4]. 
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(b) Time Histor>' of C^, and over 2 Cycles. 



Figure 5.13 Harmonic Pitching Motion about the Leading Edge 
of a 8.4% Thick Von .Mises Airfoil 
6a = 0.01 rad, o)c/Vqq = 20.0. 
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(c) Trailing Wake Pattern at the end of 2 c}'cles. 




(d) Trailing Wake Pattern at the end of 4 cycles. 



Figure 5.13 . (cont'd.) 
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(a) Pressure Distribution at a = -0.2206 rad. 



Figure 5.14 Harmonic Pitching Motion about a Pivot 0.5c ahead of 
the Leading Edge of a 8.4% Thick Von Mises Airfoil 
6a = 0.3973 rad, coc/Vqq = 0.8. 
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(b) Pressure Distribution at a = —0.0775 rad. 



Figure 5.14 . (cont'd.) 
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(c) Pressure Distribution at a = 0.1876 rad. 



Figure 5.14 . (cont'd.) 
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E. SHARP EDGE GUST FIELD PENETRATION 

1. Case- Run Definitions 

The case of an airfoil penetrating a sharp edge gust field can be computed 
using the U2DIIF code by specifying the components of gust velocity along and 
perpendicular to the freestream and the angle of attack of the airfoil. The results 
that are presented here consider the case of airfoil penetrating a sharp edge vertical 
gust at zero AOA, in view of generating information on the time-dependent lift 
resembling the classical results of Kussner [Ref 10] based on linearised theory. The 
gust front is taken to be at the leading edge at tg with only the transverse (vertical) 
component of 0.25 V qq. 

2. Results and Discussions 

Figures 5.15 and 5.16 demonstrate the variation of pressure distributions and 
trailing wake patterns respectively during and shortly after the gust front moves past 
the airfoil. It is interesting to see that the resulting wake pattern after the entire airfoil 
is submerged in the gust field is as if being split by the gust front into two portions 
curling in opposite directions. [Ref. 3] predicted a simiiar behaviour for the case of an 
undeformed gust front. Due to the use of the modified flow tangency condition to 
handle the situation when the gust front happens to fall in between tw’o nodal points, 
the pressure distributions, predicted by U2DIIF code, lie in between the results of the 
undeformed and deformed gust front models used in [Ref 3). We therefore conclude 
that this modified flow tangency condition produces sufficiently accurate results 
without adding complication in deformed gust field modeling which in turns limits the 
application to only sharp edge gusts. The present method would therefore preserve the 
generality for extending calculations to other types of gust front. Another comparison 
is made, as shown in Figure 5.17, by plotting the build-up of lift as a function of 
distance traveled by the airfoil in chord lengths. Shown in the same Figure are the 
Kussner Function and the results obtained from another case-run using a 25,5% thick 
Joukowski airfoil. 
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(a) tV'oo/c = 0.25 



Figure 5.15 Pressure Distributions at Various Time Instances Resulting 
from a 8.4% Thick Von Mises Airfoil Penetrating a 
Vertical Sharp Edge Gust of 0.25 Voq. 
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(b) tV^o/c = 0.5 



Figure 5.15 . (cont'd.) 
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Figure 5.15 . (cont'd.) 
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(d) tV'oo/c = 2.0 



Figure 5.15 . (cont’d.) 
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(e) tVoo/c = CO 



Figure 5.15 .(corn'd.) 
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Figure 5.16 Trailing Wake Patterns at Various Time Instances Resulting 
from a S.4% Thick \'on .Mises Airfoil Penetrating a 
Vertical Sharp Fdge Gust of 0.25V 
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Figure 5.16 . (cont'd.) 



101 



Kussner Function 

8.4% Thick Von Mises Airfoil 

■ X « > > ■ 25.5% Thick Joukowski .‘\irfcil 



O 




O ■■■*■■■■■ i ■ ■ . Ii ■ V m \ i ) 

0 0.5 1.0 1.5 2.0 

tV oc./c 



Figure 5.17 Normalised Time-Dependent Lift Cf. due to 
.•Xirroils of Various Thicknesses Penetrating^ 
a Vertical Sharp Edge Gust of 0.25V^. 
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V!. CONCLUDING REMARKS 



A. GENERAL COMMENTS 

The U2DIIF computer code has been developed for the purpose of 
demonstrating the successful extension of the well known panel methods, which have 
been used extensively for steady flow problems, into a powerful numerical tool for 
solving the unsteady flow problems. The mathematical modeling of the various types 
of unsteady flows has been done with the goal of preserving the generality of the 
methods. The intention is to present a method that has the minimum inherent 
limitations and restrictions so that its usage for future applications and developments 
could remain appealing. 

The validation of the U2DIIF code has been done through the various case-runs 
of the numerous types of unsteady flow problems. The results of each case-run has 
been shown to be well correlated to the results obtainable from the literature in the 
form of theoretical analyses, numerical calculations based on different variants of panel 
singularities and in cases where limited experimental data are available. The ability of 
those case-runs using an airfoil as thin as 1% thickness to produce results that 
correlate accurately to the theoretical flat plate results is perhaps a remarkable 
robustness possessed by the present unsteady flow solution methods. 

B. ENHANCING U2DIIF PROGRAM'S CAPABILITY 

It has been noted in Chapter IV that the current U2DIIF code limits the total 
number of panels to 200 and the total computation time steps to 200 also. These limits 
are not at all rigid and can be easily increased if the computer storage space is not 
critical. A point to note is that the computation time will grow rapidly as these limits 
are raised. The current linear system solver, the Gaussian elimination algorithm, used 
in the code must be concurrently improved upon to efficiently reduce the computation 
time required for the iterations in each time step. A dose examination of the matrix 
Equation 3.15, where the linear system solver is needed in every iteration, reveals that 
the coefficients of the left-hand-side matrix [ A ] are time-independent constants. 
Therefore the Gaussian elimination algorithm need only be done once for the entire 
unsteady flow calculations as far as the left-hand-side matrix is concerned. One could 
then perform, for each iteration, the manipulation of the two right-hand-sides 
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according to the steps taken for the reduction of the left-hand-side. This should cut 
down the computation time significantiy against the current method of manipulating 
both the left- and right-hand-sides simultaneously in each iteration. The savings in 
computation time w’as not pursued in the development of U2DIIF code because the 
prime concern was to demonstrate that the basic iterative solution scheme works for 
the unsteady Tow' problems. 

Another improvement that one may consider is to enable the code to be 
continuable from a time step where previous calculations w'ere terminated. One sees 
this requirement necessary not only in the case of premature termination of 
computation due to some unforeseen circumstances, but also if one needs to prolong 
the computation time. Certainly with the current code structure, one has no choice but 
to perform the calculation right from the beginning. 

.A farther extension of U2DIIF code to the computation of unsteady Row 
problems involving multiple airfoils may be w^orth pursuing. Other research w'orks that 
could be done based on U2DIIF code are in the area of incorporating more variety of 
rigid body motions into the code either in the form of closed-form equations or 
tabulated time history' of the positions and rates of motions. It is important that one 
should use as close as possible, in the code, a rigid body motion that describes the 
physical motion before generating any numerical unsteady flow results for meaningful 
comparison to test data. This fact has been well illustrated and emphasized in the 
comparison of results of case-runs involving a step change to that of a ramp change in 
AOA with a fast rate w'hich one could regard as a step in reality. However, remarkable 
difTerence in transient aerodynamics has been shown. 
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APPENDIX A 

U2DIIF PROGRAM LISTINGS 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



c 

C PROGRAM U2DIIF 

c 

C UNSTEADY MOTION OF A TWO-DIMENSIONAL AIRFOIL 

C IN INCOMPRESSIBLE INVISCID FLOW 

C USING PANEL METHODS BASED ON THE HESS & SMITH 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON /BOD/ IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+ COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON /WAR/ VYW , VXW , WAKE , DT 
COMMON /WAK2/ VYWK,VXWK 
COMMON /SING/ Q(200 ), GAMMA ,QK(200 ) ,GAMK 

COMMON /CORV/ CV( 200 ) ,XC (200 ) , YC (200 ) ,M , TD , CWX(200 ) , C'7VY(200) 
COMMON /POT/ PHI(200) ,PHIK(200) 

COMMON /GUST/ UG(200 ) , VG(200 ) , XGF , UGUST , VGUST 
COMMON /EXTV/ UE(200) 

DIMENSION XXC(200) ,YYC(200) 

C 

C INPUT FROM FILE CODE 5 AND SET UP PANEL NODES AND SLOPES 

C 



C 

C 

C 



PI = 3.1415926535 

WRITE (6,1003) 

1003 FORMAT (/////,' DATA READ FROM FILE CODE 1',//) 

CALL INDATA 
CALL SETUP 

READ (1,501) ALPI,DALP,TCON, FREQ, PIVOT, UGUST,VGUST 
WRITE (6,501) ALPI, DALP,TCON, FREQ, PIVOT, UGUST,VGUST 
501 FORMAT (7F10.6) 

READ (1,501) DELHX,DELHY, PHASE 
WRITE (6,501) DELHX,DELHY, PHASE 
READ (1,501) TF,DTS,TOL,TADJ 
WRITE (6,501) TF/DTS,TOL,TADJ 
IF (IFLAG ,EQ. 0) WRITE (6,1005) 

1005 FORMAT (///,' COORDINATES OF AIRFOIL NODES', 

+ //,3X, ' X/C ,6X, ' Y/C ,/) 

IF (IFLAG .EQ. 0) WRITE (6,1010) (X(I) ,Y(I) ,1=1 ,NODTOT+l) 
1010 FORMAT (F10.6,F10.6) 

WRITE (6,1020) SS 

1020 FORMAT(//,' AIRFOIL PERIMETER LENGTH = ',F10.6,/) 

STEADY FLOW CALCULATION AT ALPI 



ALPHA = ALPI 
WRITE (6,1030) ALPHA 

1030 FORMAT (//,' STEADY FLOW SOLUTION AT ALPHA = ',F10.6,/) 
IF (ALPHA .GT. 90.) GO TO 200 
COSALF = C0S(ilL?HA’^?I/130. ) 

3INALF = 3IM(AL?HA^PI/130. ) 

CiiLL COFISH (SINALF, COSALF) 

CALL GAUSS (1) 

CALL VELDIS (SINALF, COSALF) 

CALL FANDM( SINALF, COSALF) 

C 

C INITIALISATION FOR UNSTEADY FLOW CALCULATION TO BEGIN 



HX 


= 0.0 


HY 


= 0.0 


HXO 


= 0.0 


HYO 


= 0.0 



105 



o o o o o n o o 



nnnnot“i non non 



100 



DHX 


= 0.0 


DHY 


= 0.0 


UX 


= 0.0 


UY 


= 0.0 


ALP = ALP I 


DA 


= 0.0 


COSDA 


= 1.0 


SINDA 


= 0.0 


OMEGA 


= 0.0 


XGF 


= 0.0 


ANGLE 


= ALPI’^PI/ISO. 


COSA.NG 


= COS (.ANGLE) 


SINANG 


= SIN(ANGLE) 


DO 100 


IG = l,NODTOT 


UG(IG) 


= 0.0 


VG(IG) 


= 0.0 


PHA 


= PKASE*PI/180. 


VXW 


= COSALF 


VYW 


= SINALF 


GAMK 


= GAMMA 


T 


= 0.0 


M 


= 0.0 


TOLD 


= 0.0 



+ ATAN(VGUST/(1 .+UGUST) ) 



RIGID BODY MOTIONS OF AIRFOIL 

IF (FREQ .NE. 0.0) GO TO 1 

IF (DALP .EO. 0.0) GO TO 2 

IF (TOON .Ng. O.O) GO TO 3 

ALPHA = ALP I + DALP 
COSALF = COS (ALPHA ^^'PI/ 130. ) 

SINALF = SIN(ALPHA’^PI/180. ) 

3 DT = DTS 

TD = DTS 

GO TO 60 

2 IF ((UGUST .EQ. 0.0) .AND. (VGUST .EQ. 0.0)) GO TO 200 
DT = DTS 

TD = DTS 

GO TO 60 

1 DT = 2.0*PI/(FREQ*DTS) 

TD = DT 

60 T = DT 

WRITE (6,1051) 

10^1 FORMAT (///// ' :k:k:k:k:k:k:k:k'k'k'k'k'k'k'k:k:k'k:k:k:k:k:k:k:k:k:k:k:k:k'k:k:k:k:k:k:k:kii I / 

+ I ***' begin unsteady flow solution ****'/, 

-I- I 'k:k:k:k:k'k7^:k'k:k:ki(:i(:i>c:kii:k:k:k:k:k7K:ki^:k:k:k:k:k'k'k'k:k'k:k"k'k:k:k I J 

40 M = M + 1 

IF (T .GT. TF) GO TO 200 

STORE CORE VORTEX COORDINATES FOR TIME STEP ADJUSTMENTS 

IF (M .EQ. 1) GO TO 50 
DO 51 I = 1,M-1 

xxcU) = XC(I) 

51 YYC(I) = YC(I) 

50 IF (FREQ .NE. 0.0) GO TO 11 
IF (DALP .EO. 0.0) GO TO 22 
IF (TCON .me. 0.0) GO TO 33 

STEP CHANGE IN AOA 

IF (TADJ .NE. 0.0) GO TO 70 
TD = FLOAT (M+l)*DTS 

GO TO 70 



MODIFIED RAMP CHANGE IN AOA 

33 IF (T .GT. TCON) GO TO 34 

DAL = DALP * (3.-2.*T/TCON)*(T/TCON)**2 

ALPHA = ALP I + DAL 
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on non non 



COSALF = COS(ALPHA^PI/180. ) 

SIMALF = SIN(AL?HA’^?I/130. ) 

DA = ALPHA - ALP 

COSDA = COS(DA*PI/130. ) 

SINDA = SIM(DA’^PI/130. ) 

OMEGA = - (DALP’^PI/130. ) (5 . ’^T/ (TCON’^TCON) ) * (1,-T/TCON) 

DHX = PIVOT (1.- COSDA) 

DKY = - PIVOT ^ SINDA 

UY = PIVOT ^ OMEGA 

MTCOM = M 

GO TO 70 

34 DAL =0.0 

. ALPHA = ALP I + DALP 

COSALF = COS(ALPHA*PI/180.) 

SINALF = SIN(ALPHA*PI/180. ) 

DA =0.0 

COSDA =1.0 

SINDA =0.0 

OMEGA =0.0 

DHX =0.0 

DHY =0.0 

UY =0.0 

IF (TADJ .ME. 0.0) GO TO 70 
TD = FLOAT (M+l-MTCON)’^DTS 

GO TO 70 

SHARP EDGE GUST (UGUST AND/OR VGUST) 

22 XGF = T 

DO 110 IG = l,NODTOT 

ug(:g) =0.0 

VG(ZG) =0.0 

XG = X(IG)^COSALF + Y ( IG) ’^SINALF 

XGPl = X(IG+1) ^COSALF + Y( IG+1 ) =*^SINALF 

IF (IG .LT. NLOWER+1) GO TO 120 
IF (XGF .LE. XG) GO TO 110 
IF (XGF .GE. XGPl) GO TO 111 
FAC = (XGF - XG)/(XGP1 - XG) 

UG(IG) = UGUST’^FAC 
VG(IG) = VGUST^FAC 
GO TO 110 

111 UG(IG) = UGUST 
VG(IG) = VGUST 

GO TO 110 

120 IF (XGF ,LE. XGPl) GO TO 110 

IF (XGF .GE. XG) GO TO 121 

FAC = (XGF - XGP1)/(XG - XGPl) 

UG(IG) = UGUST’^FAC 
VG(IG) = VGUST’^FAC 
GO TO 110 

121 UG(IG) = UGUST 

VG(IG) = VGUST 

110 CONTINUE 

IF (XGF .LE. COSALF) MGUST = M 
IF (TADJ .NE. 0.0) GO TO 70 

IF (XGF .GT. COSALF) TD = FLOAT (M+1 -MGUST )*DTS 
GO TO 70 

TRANSLATION HARMONIC OSCILLATION 

11 IF (DALP .NE. 0.0) GO TO 12 

HX = DELHX ^ SIN(FREQ*T + PHA) 

HY = DELHY * SIN(FREQ*T) 

DHX = HX - HXO 

DHY = HY - HYO 

UX = DELHX=*^FREQ*COS(FREQ*T+PHA) 

UY = delhy*freq*cos(freq*t) 

GO TO 70 

ROTATIONAL HARMONIC OSCILLATION 
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12 



C 

C 

c 



70 



90 

80 

1001 

1004 



C 

c 



DAL = DAL?’^SIN(?RE0^T) 

ALPHA = ALPI T DAL 

COSALF = COS(ALPKA’^?I/1SO. ) 

SINAL? = SIN(AL?KA^PI/180. ) 

DA = ALPHA - ALP 

COSDA = COS(DA’^PI/130.) 

SIMDA = SIN(DA^PI/130. ) 

OMEGA = - ('DALP^PI/180. ) * FREQ COS(FREQ*T) 

UY = PIVOT * OMEGA 

DHX = PIVOT * (1. -COSDA) 

DHY = - PIVOT ^ SIMDA 



TRAMSFORM CORE VORTEX COORDINATES W. R. T, NEW AIRFOIL POSITION 



IF (M .EO. 1) GO TO 30 

DO 90 'I = 1,M-1 

XC(I) = XXC(I) + CWX(I) ^ DT 

YC(I) = YYC(I) + C'/'7Y(I) ^ DT 

XCO = XC ( I ) 

YCO = YC(I) 

XC(I) = XCO’^COSDA - YCO^SINDA + DHX 

YC(I) = XCO^SINDA + YCO^COSDA + DHY 
CONTINUE 

WRITE (6,1001) T,DT 

FORMAT llllli,' TIME STEP TK = ‘ , FiO . 6 , lOX , ' TK - 
WRITE (6,1004) ALPHA , OMEGA , UX , UY 
FORMAT (/,' ALPHA(T) = ',F10.6,5X,' OMEGA(T) = 

+ ' U(T) = ' ,F10.6,5X, ' V(T) = ' 

+ IX,' MITR VXW VYW WAKE THETA 



TKMl = ' ,F10.6,/) 



FIO, 6,/, 
,F10.6,///, ^ 

GAMMA' ,/) 



CALCULATE THE TRAILING EDGE WAKE ELEMENT 



MITR = 0 

10 'WAKE = SQRT('-7YW^VYV7+VXV?’^VXW)=^DT 
THENP 1 = ATAN2 ( VYW , VXW ) 

COSTHE(NPl) = COS(THENPl) 

SINTHE(MPl) = S IN ( THENP 1) 

WRITE (6,1002) NITR, VXW, VYW, WAKE, THENP 1,GAMK 
1002 FORIIAT (I5,4F10.6,E14.6) 

X(NP2) = X(NPl) + WAKE’^COSTHE(NPl) 

Y(mP 2) = Y(NPl) + WAKE’''SINTHE(NP1) 

CALL INFL (NITR) 

CALL COEF (SINALF, COSALF, OMEGA, UX,UY) 

CALL GAUSS (2) 

CALL KUTTA (ALPHA, SINALF, COSALF, OMEGA, UX,UY) 

CALL TEWAK (SINALF , COSALF) 

TOLl = ABS(VYW - VYWK) 

TOL2 = ABS(VXW - VX^-JK) 

IF ((TOLl ,LT. TOL) .AND. (T0L2 .LT. TOL)) GO TO 20 

VYW = VYWK 

VXW = VXWK 

IF (NITR .GT. 200) STOP 

NITR = NITR + 1 

GO TO 10 

20 WRITE (6,1011) NITR 

1011 FORMAT < ! I CONVERGED SOLUTION OBTAINED AFTER NITR = ',12) 

:ALL PRESS (SINALF, COSALF, OMEGA, UX,UY) 

IF ' (UGUST .EO. 0.0) .AND. (VGUST . EQ , 0.0)) GO TO 300 
CALL FANDM (SiNANG , COSAI'IG) 

GO TO 400 

300 CALL FANDM (SINALF , COSALF) 

400 CONTINUE 
C 

C ADJUST TIME STEP (TADJ .NE. 0.0) IF NECESSARY 
C 

IF (TADJ .EQ. 0.0) GO TO 95 
WRITE (5,2001) 

2001 FORMAT (//,' DO YOU WANT TO ADJUST TIME STEP ? 0 - NO, 1 - YES') 
READ (5,*) IDT 
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IF (IDT .EQ. 0) GO TO 95 
DT = TADJ * DT 

T = TOLD + DT 

WRITE (6,1006) 

1006 FORMAT (//,' BACK-TRACK COMPUTATION AND ADJUST TIME-STEP',//) 
GO TO 50 



C 

C 

c 



WAKE ELEMENT LEAVES TRAILING EDGE AS A CORE-VORTEX 



95 



CV(M) 

XC(M) 

YC(M) 



= SS *( GAMMA -GAMK) 



= X(NPl) 
= Y(NPl) 



CWX(M) = VXW 
CWY(M) = VYW 
WRITE " 



+ 0.5*WAKE’^COSTHE(NP1) 
+ 0.5’^WAKE’^SINTHE(NP1) 



'(6,1052) 

1052 FORMAT (//,' TRAILING VORTICES DATA',//, 
+ 4X, 'M' ,4X, 'X(M) ' ,6X, 'Y(M) ' ,6X, 'CIRC ,/) 
DO 900 I = 1,M 

WRITE (6,1050) I,XC(I) ,YC(I) ,CV(I) 
FORMAT(I5,3F10.6) 

CALL GORVOR (SINALF , COSALF) 



900 

1050 



C 

C 

C 



30 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 



RE-INITIALISE PARAMETERS FOR NEXT TIME STEP CALCULATION 



DO 30 


I 


= l,NOD 


Q(I) , 


= 


QK(I) ^ 


PHI ( I ) 




PHIK(I) 


CONTINUE 






GAMMA 


= 


GAMK 


ALP 


= 


ALPHA 


HXO 




HX 


HYO 




HY 


TOLD 


— 


T 


DT 




TD 


T 




T + TD 


GO TO 40 





200 STOP 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



SUBROUTINE BODY(Z,SIGN,X,Y) 

RETURN COORDINATES OF POINT ON THE BODY SURFACE 

Z = NODE-SPACING PARAMETER 
X,Y = CARTESIAN COORDINATES 
SIGN = +1. FOR UPPER SURFACE 
-1. FOR LOWER SURFACE 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE BODY(Z,SIGN,X,Y) 

COMMON /PAR/ NACA,TAU,EPSMAX,PTMAX 
IF (SIGN .LT. 0.0) Z = 1. - Z 
CALL NACA45 ( Z , THICK , CAMBER , BETA ) 

X = Z - SIGN*THICK’''SIN(BETA) 

Y = CAMBER + SIGN’''THICK’''COS(BETA) 

RETURI'I 

END 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C SUBROUTINE COEF (SINALF, COSALF, OMEGA, UX,UY) 

C 

C SET COEFFICIENTS OF N EQUS ARISING FROM FLOW 

C TANGENCY CONDITIONS AT MID POINTS OF PANELS 

C SOLVING THE N- SOURCE STRENGTHS IN TERMS OF THE 

C VORTICITY STRENGTH (RESULTING IN 2 RHS) 

C KUTTA CONDITION IS SATISFIED SEPARATELY TO OBTAIN 

C THE VORTICITY STRENGTH 

C THIS SOLUTION METHOD IS DESIRED FOR UNSTEADY FLOW 

C 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE COEF (SINALF , COS AL? , OMEGA , UK , UY ) 

COMMON /BOD/ IFLAG,NLOWER,NU??ER,NODTOT,X(202) ,Y(202) , 

+ COSTHE(201) ,SINTHE(201) ,SS,MP1,NP2 

COMMON /COF/ A(201,2ll) ,NEQS 
COMMON /SING/ Q ( 200 ), GAMMA, QK( 200) ,GAMK 
COMMON /WAK/ VYW , VXW , WAKE , DT 

COMMON /CORV/ CV(200) ,XC(200) ,YC(200) ,M,TD,CCVX(200) ,CCVY(200) 

COMMON /INFl/ AAN(201,201) ,BBN(201,201) ,AY!^IP1(201),3YN?1(201) 

COMMON /INF2/ CCN( 201 , 200 ) , CCT ( 201 , 200 ) , CYNPl ( 200) , CI'^NPl (200) 

COMMON /GUST/ UG( 200 ) , VG( 200 ) , XGF , UGUST , VGUST 

NEQS = NODTOT 

NPl = NODTOT + 1 

NP2 = NODTOT + 2 



INITIALISE COEFFICIENTS 



DO 90 I = 1, NODTOT 

DO 90 J = 1,NP2 

90 A(I,J) = 0.0 

SET LHS MATRIX A(I,J) 



DO 120 I = 1, NODTOT 

Xl^ID = 0.5 * (Xh) + X(I-fl)) 

YMID = 0.5 * (Y(I) + Y(I+1)) 

B = 0.0 

DO 110 J = 1, NODTOT 

A(I,J) = AAN(I,J) 

B = B + BBN(I,J) 

110 CONTINUE 



FILL IN THE RIGHT HAND SIDE 

A(I,NP1) = -B + BBM(I,NP1)*SS/WAKE 
A(I,NP2) = -BBN{I,NP1)*GAMMA*SS/WAKE 
+ + SINTHE(I)* ((l.+UG(l))*COSALF-VG(I)*SINALF+UX) 

+ - COSTHE(I)* ((1.+UG(I))*SINALF+VG(I)*C0SALF+UY) 

+ + OMEGA*(YMID’^SINTHE(I) + XMID’^COSTHE (I ) ) 

ADD CORE VORTEX CONTRIBUTION 



IF (M .EQ. 1) GO TO 
Ml^l = M - 1 

1,MM1 



140 



DO 100 N 

A(I,NP2) = A(i,NP2) - CCN(I,N)*CV(N) 

100 CONTINUE 
140 CONTINUE 
120 CONTINUE 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



SUBROUTINE COFISH(SINALF , COSALF) 

SET COEFFICIENTS OF LINEAR SYSTEM - M+1 EQUATIONS 

M EOUS - FLOW T.ANGENCY AT MID POINTS ~0F PANELS 
1 EdU - X.ATTA CONDITION AT TRAILING EDGE PANELS 
THIS SOLUTION METHOD IS EFFECTIVE FOR STE.ADY FLOW, MO 
ITERATION IS REQUIRED, N-SOURCE STRENGTHS AND 1 
VORTICITY STRENGTH ARE SOLVED SIMULTANEOUSLY 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE COFISH(SIN.ALF , COSALF) 

COMMON /BOD/ IFLAG , NLOWER , NUPPER ,NODTOT , X( 202) . Y( 202) , 

+ COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON /COF/ A(201,211) ,KUTTA 
COMMON /MUM/ PI,PI2INV 
KUTTA = NODTOT + 1 



no 



n n o o n o r:> n n n 



INITIALISE COEFFICIENTS 



C 

c 



90 



C 

c 

c 



c 

c 

c 



100 



c 

c 

c 

c 
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DO 90 J = 1,KUTTA 
A(KUTTA,J) =0.0 

SET VN = 0 AT MID-POINT OF I-TH PANEL 

DO 120 I = l,NODTOT 

XHID = .5-(X(I) + X(I + 1)) 

YMID = .5=^(Y(I) + Y(I+I)) 

A(I,XUTTA) =0.0 

FIND CONTRIBUTION OF J-TH PANEL 

DO 110 J = l^NODTOT 
FLOG =0.0 

FTAN = PI 

I? (J .EO. I) GO TO 100 

DKJ = XMID - X(J) 

DXJ? = XMID - X(J+1) 

DYJ = YMID - Y(J) 

DYJP = YMID - Y(J+1) 

FLOG = .5^ALOG( (DXJ?’^DXJP+DYJP^DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
FTAN = ATAN2(DYJP^DXJ-DXJP^DYJ,DXJP^DXJ+DYJP^DYJ) 

CTIMTJ = COSTHE( I)’^COSTHE( J) + 5INTHE ( I ) ^SINTHE ( J ) 

5TIMTJ = SINTKE(I)*COSTHE( J) - COSTHE ( I ) ^SINTHE ( J) 

A(I,J) = PI2INV’^(FTAN^CTIMTJ + FLOG^STIMTJ) 

3 = ?I2INV*(FL0G*CTIMTJ - FTAN’^STIMTJ) 

A(I,XUTTA) = A(I,KUTTA) + B 

IF ((I .GT. 1) .AND. (I .LT. NODTOT))GO TO 110 

I? I-TH PANEL TOUCHES TRAILING EDGE, 

ADD CONTRIBUTION TO XUTTA CONDITION 

A(KUTTA,J) = A(KUTTA,J) - 3 
A(KUTTA,KUTTA) = A(KUTTA,KUTTA) + A(l,J) 

CONTINUE 



C 

C FILL IN KNOWN SIDES 

C 

A(I,KUTTA+1) = SINTHE(I)*COSALF - COSTHE(I )*SINALF 
120 CONTINUE 

A(KUTTA,KUTTA+1) = - (COSTHE(l) + COSTHE (NODTOT)) *COSALF 
+ - (SINTHE(l) + SINTHECNODTOTH’^SINALF 

RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C SUBROUTINE CORVOR (SINALF , COSALF) 

C 

C COMPUTE THE LOCAL VELOCITIES OF CORE VORTICES 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE CORVOR (SINALF , COSALF) 

COMMON /BOD/ IFLAG , NLOWER , NUPPER ,NODTOT , X(202 ) , Y(202 ) , 

+ COSTHE(201) ,SINTKE(201) ,SS,NP1,NP2 

COMMON /SING/ 0 ( 200 ), GAMMA , OK ( 200 ), GAMK 
COMMON /WAK/ 7YW , VXW , WAKE , DT 

COMMON /CORV/ CV(200) ,XC(200) ,YC(200) , M ,TD , CCVX( 200 ) , CC7Y( 200 ) 
COMMON /INF3/ AMY( 200 , 201 ) , BMY( 200 , 201 ) 

COMMON /INF4/ CMY(200 , 200) , CMX(200 , 200) 

COMMON /POT/ PHI(200) ,PHIK(200) 

COMMON /GUST/ UG(200 ) , VG( 200 ) , XGF , UGUST , VGUST 
IF (M .EQ. 1) GO TO 40 
MMl = M - 1 



C 

C VELOCITY COMPONENTS OF CORE VORTICES AT CURRENT TIME STEP 

C 



UGC 

VGC 



0.0 

0.0 
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c 

c 

c 



20 



DO 10 
XG 

I? (XG 
UGC 
VGC 
VY 

-S- 

’ vx 

DO 20 

'Pi 

7X 

CONTINUE 



N = 1,MM1 

= XC(N)^COSALF + YC(M)*SINALF 
,GT, XGF) GO TO 5 
= UGUST 
= VGUST 

= S S 3MY ( N , N? 1 ) ’^ ( GAMMA - G AMK ) / WAKE + 

( 1 , +UGC) ^SINALF+VGC’^COSALF 
= SS*A^rf (N , NP 1 ) ^ ( GAMMA-GAMX ) /WAKE+ 

(1 .^UGC)^COSALF-VGC’^SINAL 
J = I,NODTOT 

= VY + AMY(N, J)’^OK(J) + BMY ( N , J ) ’^GAMK 
= VX - BMY(N, J)^6K(J) + AMY(N, J)’^GAMK 



ADD CORE VORTEX CONTRIBUTION 



30 



C 

C 

C 



DO 
Ir 
VY 
VX 

CONTINUE 



30 MC = 1,MM1 
(MC .EO. N) GO TO 30 



= 'VY + CMY(N,MC)’^CV(MC) 
= VX + CMX(N,MC)*CV(MC) 



COORDINATES 

CCVX(N) = VX 
^7Y 



OF CORE VORTICES AT NEXT TIME STEP 





CC"7Y(N) 


10 


CONTINUE 


40 


CONTINUE 


50 


CONTINUE 




RETURN 




END 


^ 


W v.« ^ ^ W Vw ^ ' 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c c 

C SUBROUTINE FANDM(SIMALF,COSALF) C 

C C 

C COMPUTE AND PRINT OUT CD, CL, CM C 

C INTEGRATE PRESSURE DISTRIBUTION BY TRAPEZOIDAL RULE C 

C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE FANDM(SINALF , COSALF) 

COMMON /BCD/ IFLAG ,NLOWER , NUPPER ,NODTOT , X(202 ) , Y(202 ) , 

+ COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON /CPD/ CP (200) 



100 



CFX 
CFY 
CM 

DO 100 

XMID 

YMID 

DX 

DY 

CFX 

CFY 

CM 

CONTINUE 

CD 



0.0 
0.0 
0.0 

= l,NODTOT 
.5’^(X(I) + X(I+1)) 

.5’^(Y(I) + Y I+l)) 

X(I+1 - XI 

Y(l+l) - Y(I) 

CFX + CP(I)’^DY 
CFY - CP(I)^DX 

CM + CP(I)^(DX*XMID + DY*YMID) 



CM =' ,F10.6) 



CFX’^ COSALF + CFY*SIN.ALF 
CFY^COSALF - CFX=^SINALF 
WRITE <'6,1000) CD, CL, CM 

1000 FORMAT(//,' CD =',F10.o,' CL =',F10.6,’ 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

C SUBROUTINE GAUSS (NRHS) C 

c c 

C SOLUTION OF LINEAR ALGEBRAIC SYSTEM BY C 

C GAUSS ELIMINATION WITH PARTIAL PIVOTING C 

c c 

C A = COEFFICIENT MATRIX C 

C NEQNS = NUMBER OF EQUATIONS C 
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C NRHS = NUMBER OF RIGHT HAND SIDES 

C 

C RIGHT-HAND SIDES AND SOLUTIONS STORED IN 

C COLUMNS NEQNS+1 THRU NEQNS+NRHS OF A 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE GAUSS (NRHS) 

COMMON /COF/ A(201,211) ,NEQNS 
NP = NEONS + 1 

NTOT = NEONS + NRHS 

C 

C GAUSS REDUCTION 

C 

DO 150 I = 2,NEQNS 
C 

C — SEARCH FOR LARGEST ENTRY IN (I-l)TH COLUMN 

C ON OR BELOW MAIN DIAGONAL 

C 

IM =1-1 

IMAX = IM 

AMAX = ABS(A(IM,IM)) 

DO 110 J = I, NEONS 

IF (AMAX .GE. AB§(A(J,IM))) GO TO 110 
IMAX = J 

AMAX = ABS(A(J,IM)) 

110 CONTINUE 
C 

C — SWITCH (l-l)TH AND IMAXTH EQUATIONS 

C 

IF (IMAX .EQ. IM) GO TO 140 
DO 130 J = IM,NTOT 
TEMP = A(IM,J) 

A(IM,J) = A(IMAX,J) 

A(IMAX,J) = TEMP 
130 CONTINUE 



C 

C ELIMINATE (I-l)TH UNKNOWN FROM 

C ITH THRU (NEQNS)TH EQUATIONS 

C 

140 DO 150 J = I,NEQNS 

R = A(J,IM)/A(IM,IM) 

DO 150 K = I, NTOT 

150 A(J,K) = A(J,K) - R^A(IM,K) 

C 

C BACK SUBSTITUTION 

C 



200 

210 

220 



DO 220 K = NP,NTOT 

A(NEQNS,K) = A(NEQNS,K)/A(NEQNS,NEQNS) 
DO 210 L = 2,NEQNS 
I = NEQNS + 1 - L 

IP =1+1 

DO 200 J = IP, NEQNS 
k(l,K) = A(I,K) - A(I J)’*'A(J,K) 

A(I,K) = AU,K)/A(I,I) 

CONTINUE 

RETURN 

END 



SUBROUTINE INDATA 



C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



SET PARAMETERS OF BODY SHAPE 

FLOW SITUATION, AND NODE DISTRIBUTION 

USER MUST INPUT 

NLOWER = NUMBER OF NODES ON LOWER SURFACE 
NUPPER = NUMBER OF NODES ON UPPER SURFACE 
PLUS DATA ON BODY AND SUBROUTINE BODY 
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nnnnnnn nnnnnn 



SUBROUTINE INDATA 
DIMENSION TITLE (20) 

COMMON /BOD/ IFLAG ,NL0WER ,NUPPER ,N0DT0T , X( 202 ) , Y( 202) , 

+ COSTHE(201) ,SINTKE(201) ,SS,NP1,NP2 

COMMON /PAR/ NACA , TAU , EPSMAX , PTMAX 
READ (1,501) ITITLE 
WRITE (6,501) ITITLE 
DO 10 I = 1, ITITLE 
READ (1,502) TITLE 

10 WRITE (6,503) TITLE 

501 FORMAT (31 5) 

502 FORMAT (20A4) 

503 FORMAT (IX, 20A4) 

READ (1,501) IFLAG,NLOWER,NUPPER 

WRITE (6,501) IFLAG,NLOWER,NUPPER 

IF UfLAG .NE. 0) RSTUPJ'I 

READ (1,501) NACA 

WRITE (6,501) NACA 

lEPS = NACA/ 1000 

IPTMAX = NACA/100 - 10*IEPS 

ITAU = NACA - lOOO’^IEPS - lOO’^IPTMAX 

EPSMAX = lEPS^O.Ol 

PTMAX = IPTMAX’^O.l 

TAU = ITAU^O.Ol 

IF (lEPS .LT. 10) RETURN 

PTl^X = 0.2025 

EPSMAX = 2.6595=^PTMAX**3 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE INFL (MITR) 

CALCULATE INFLUENCE COEFFICIENTS 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE INFL (NITR) 

COMMON /BOD/ IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+ COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON /NUM/ PI,PI2INV 
COMMON /WAK/ VYW , VXV7, WAKE , DT 

COMMON /CORV/ CV(200) ,XC(200) ,YC(200) ,M,TD,CCVX(200) ,CCVY(200) 
COMMON /INFl/ AAN(201,201) ,BBN(201,201) ,AYNP1(201) ,BYNP1(201) 

“ ■ ‘ ,CCT(201,200) ,CYNP1(200) ,CXNP1(200) 

,BMY(200,201) 

,CMX(200,200) 

.GT. 0)) GO TO 510 



COMMON /INF2/ CCN(201,200' 
COMMON /INF3/ AMY( 200,201' 
COMMON /INF4/ CMY(200,200' 
IF ((M .GT. 1) .OR. (NITR' 



AAN(I,J) : NORMAL VELOCITY INDUCED AT MID-POINT OF I-TH PANEL 

BY UNIT STRENGTH DISTRIBUTED SOURCE ON THE J-TH PANEL 

BBN(I,J) : NORMAL VELOCITY INDUCED AT MID-POINT OF I-TH PANEL 

BY UNIT STRENGTH DISTRIBUTED VORTEX ON THE J-TH PANEL 



100 



= l,NODTOT 
.5*(X(I) + X(I+1)) 
.5’='(Y(I) ^ Y(I + 1)) 



= 0.0 
= PI 
.EQ. I) 

= XMID 
= XMID 
= YMID 
= YMID 



MODTOT 



DO 120 
XMID 
YMID 
DO 110 
FLOG 
FTAN 
IF (J 
DXJ 
DXJP 
DYJ 
DYJP 
FLOG 
FTAlvI 
CTIMTJ 
STIMTJ 

AAI'I(I,J) = Pl2lNV*(FTAN’^CTiMTJ + FLOG^STIMTJ) 



GO TO 100 
X(J) 

X(J+1) 

Y J) ^ 

Y( J+1) 



5^ALOG( (DXJP*DXJP+DYJP’^DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
= ATAN2(DYJP*DXJ-DXJP’^DYJ,DXJP*DXJ+DYJP’^DYJ) 

= COSTHE(i)’^COSTHE(J) + SINTHE(I)^SINTHE(J) 

= SINTHE(I)’^COSTHE(J) - COSTHE ( I ) *SINTHE ( J) 
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110 

120 

510 



135 



BBN(I,J) = ?I2INV^(FL0G*CTIMTJ - FTAN^STIMTJ) 

CONTINUE 
CONTINUE 
CONTINUE 
I = NPl 

XMID = .5=*'(X(I) + X(I+i)) 

YMID = .S^{Y(I} + T(I+1)) 

DO 130 J = 1,NP1 

FLOG =0.0 

FTAN = ?I 

IF (J .SO. I) GO TO 135 
DXJ = Xl^ID - X(J) 

DXJP = XMID - X(J+1) 

DYJ = YMID - Y(J) 

DYJP = YMID - Y(J+1) 

FLOG = .5*aLOG((DXJP’^DXJP+DYJP’^DYJP)/(DXJ*DXJ+DYJ^DYJ)) 

FTAN = ATAN2(DYJP’^DXJ-DXJP’^DYJ,DXJP’^DXJ+DYJP=*'DYJ) 

CTIMTJ = COSTHE(I)’^COSTHE(J) + SINTHE ( I ) ’^SINTHE ( J) 

STIMTJ = 5INTHE(I)^C0STHE(J) - COSTHE ( I ) ’^SINTHE ( J) 
AAN(I,J) = PI2INV*(FTAN’^CTIMTJ + FLOG’^STIMT J) 

•BBN(I,J) = PI2INV^(FLOG^CTTMTJ - FTAN^STIMTJ) 



AYNPl(J) : Y - VELOCITY INDUCED AT MID POINT OF WAKE ELEMENT 
(NPl-TK PANEL) BY UNIT STRENGTH DISTRIBUTED SOURCE 
ON J-TH PANEL 



BYNPl(J) : Y - VELOCITY INDUCED AT MID POINT OF WAKE ELEMENT 
(NPl-TH PANEL) 3Y UNIT STRENGTH DISTRIBUTED VORTEX 
ON J-TH PANEL 



130 



140 



AYNPl(J) 
3YNP1(J) 
CONTINUE 
DO 140 
XIHD 
YMID 
J 

DXJ 
DXJP 
DYJ 
DYJP 
FLOG 
FTAN 
CTIMTJ 
STIMTJ 
AAN(I, J) 
BBN(I, J) 
CONTINUE 



= ?I2INV*(FTAN*COSTHE(J) - FLOODS INTHE (J) ) 
= PI2IMV*(FLOG^COSTHE(J) + FTAN^SINTHE ( J) ) 



I = i,NODTOT 
= .5’^(X(I) + X(I+1)) 
= .5*(Y(I) + Y(I+1)) 



= NPl 

= XMID - X(J) 

= XMID - X(J+1) 

= YMID - Y(J) 

= YMID - Y(J+1) 

= .5’^ALOG( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
= ATAN2 ( DY JP*DXJ-DXJP*DY J , DXJP^DXJ+DY JP*DYJ ) 

= COSTHE(I)*COSTHE(J) + SINTHE ( I ) ’^SINTHE ( J) 

= SINTHE(I)’^COSTHE(J) - COSTHE ( I ) ’^SINTHE ( J) 

= PI2INV*(FTAN’^CTIMTJ + FLOG^STIMTJ) 

= PI2INV*(FLOG*CTIMTJ - FTAN’^STIMTJ) 



CCN(I,J) : NORMAL VELOCITY INDUCED AT MID-POINT OF I-TH PANEL 
BY UNIT STRENGTH N-TH CORE VORTEX 

CCT(I,J) : TANGENTIAL VELOCITY INDUCED AT MID-POINT OF I-TH PANEL 
BY UNIT STRENGTH N-TH CORE VORTEX 

IF (M .EQ. 1) RETURN 

MMl = M - 1 

IF (NITR .GT. 0) GO TO 520 

DO 220 I = i.MODTOT 

XMID = 0.5’^(X(I) + X(I+1)) 

YMID = 0.5*(Y(I) + Y(I+1)) 

DO 210 N = 1,MM1 
DX = XMID - XC(N) 

DY = YMID - YC(N) 

DIST = SQRT(DX*DX+DY’^DY) 

COSTHN = DX/DIST 
SINTHN = DY/DIST 

CTIMTN = COSTHE(I)^COSTHN + SINTHE ( I ) ’^SINTHN 
STIMTN = SINTHE(i)’^COSTHN - COSTHE ( I ) ’^SINTHN 
CCN(I,N) = -PI2INV*CTIMTN/DIST 
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CCT(I,N) = -?I2INV*STIMTN/DIST 
210 CONTINUE 
220 CONTINUE 
520 CONTINUE 

: = NPl 

XHID = 0.5^(X(I) + X(I+1)) 

YMID = 0.5^(Y(I) + Y(I+1)) 

DO 230 N = 1,MM1 
DX = XHID - XC(N) 

DY .= YHID - YC(N) 

DIST = S0RT(DX^DX+DY*DY) 

COSTHN = DX/DIST 
SINTKN = DY/DIST 

CTIMTN = C0STHE(I)*C0STHN + SINTHE(I)’^SINTHN 
STIMTN = SINTKE(I)^COSTHN - COSTKE ( I ) *SINTHN 
CCN(I,N) = -?I2INV^CTIHTN/DIST 
CCT(I,N) = -PI2INV’^STIHTN/DIST 

CYNPl(N) : Y - VELOCITY INDUCED AT MID POINT OF WAKE ELEMENT 
(NPl-TH PANEL) BY UNIT STRENGTH N-TH CORE VORTEX 

CXNPl(N) : X - VELOCITY INDUCED AT MID POINT OF WAKE ELEMENT 
(NPl-TH PANEL) BY UNIT STRENGTH N-TH CORE VORTEX 

CYNPl(M) = -PIYINV’^COSTKN/DIST 
CXMPl(N) = +?I2IMV’^SINTHN/DIST 
230 CONTINUE 

AMY(N,J) : Y - VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 
STRENGTH DISTRIBUTED SOURCE ON THE J-TH PANEL 

3ITY(N,J) : Y - VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 
STRENGTH DISTRIBUTED VORTEX ON THE J-TH PANEL 



IF (NITR , 
DO 320 N 
XMID = 
YMID 

DO 310 J 
DXJ = 

DXJ? 

DYJ = 

DYJP = 
FLOG 

FTAN = 
AMY(N,J) = 
BMY(N,J) = 
310 CONTINUE 
320 CONTINUE 
530 CONTINUE 

DO 330 N 

XMID = 

YMID = 

J 

DXJ 

DXJ? 

DYJ 

DYJ? 

FLOG = 

FTAN 

AIf2(N,J) = 
BMY(M, J) 
330 CONTINUE 

C^^Y(N,MC) 



GT. 0) GO TO 530 

= 1,MM1 

XC(M) 

YC(N) 

= 1,N0DT0T 
XI^ID - X(J) 

XI^ID - X(J+1) 

YMID - Y(J) 

YMID - Y(J+1) 

. 5*ALOG ( ( DX JP*DX JP+DY JP*DY JP ) / ( DX J’^DXJ+DY J*DY J ) ) 
ATAN2 ( DYJP*DXJ-DXJP*DYJ , DXJP^DXJ+DYJP^DYJ) 

: PI2INV*(FTAN*C0STHE(J) - FLOG’^SINTHE ( J^ 

: PI2INV*(FL0G*C0STHE(J) + FTAN^SINTHE ( J) ) 



= 1,MM1 
: XC(N) 

' YC(N) 

: NPl 

= XMID - X(J) 

: XMID - X(J+1) 

: YMID - Y(J) 

• YMID - Y(J-i-l) 

: . 5'AL0G( (DXJP’^DXJP^DYJP^DYJ?)/ (DXJ^DXJ+DYJ^DYJ) ) 
: ATAN2 ( DY JP’^DX J -DX JP^DY J , DX JP’^DX J+DY JP^DYJ ) 

= PI2INV*(FTAM*COSTHE(J) - FLOG’^SINTHE ( J) ) 

= PI2INV^(FL0G*C0STHE(J) + FTAN^SINTHE ( J) ) 



CMX(N,MC) 



: Y - VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 
STRENGTH MC-TH CORE VORTEX OTHER THAN ITSELF 

: X - VELOCITY INDUCED AT N-TH CORE VORTEX BY UNIT 
STRENGTH MC-TH CORE VORTEX OTHER THAN ITSELF 
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IF (NITR .GT. 0) RETURN 
DO 420 N = 1,MM1 
XMID = XC(N) 

YMID = YC(N) 

DO 410 MC = 1,MM1 
IF (NC .EQ. N) GO TO 410 
DX = XMID - XC(MC) 

DY = YMID - YC(MC) 

DIET = SORT(DX*DX+DY’^DY) 

COSTHN = DX/DIST 
SINTHN = DY/DIST 
CMY(N,MC) = -PI2INV’^C0STHN/DIST 
CMX(N,MC) = +PI2INV*SINTHN/DIST 
410 CONTINUE 
420 CONTINUE 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE KUTTA ( ALPHA, SINALF, COSALF, OMEGA, UX,tTY) 

USING KUTTA CONDITION TO DETERMINE VORTICITY 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



SUBROUTINE KUTTA (ALPHA , SINALF , COSALF , OMEGA , UX , UY) 

COMMON /BOD/ IFLAG,NLOWER,NUPPER,NODTOT ,X(202) ,Y(202) , 
COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 
COMMON /COF/ A(201,211) ,NEQS 
COMMON /SING/ Q ( 200 ), GAMMA, QK( 200) ,GAMK 
COMON /WAK/ VYW,VXW,WAKE,DT 

COMMON /CORV/ CV(200) ,XC(200) ,YC(200) , M, TD , CCVX ( 200 ) , CCVY( 200 ) 
COMI'ION /INFl/ AAN(201 ,201) ,BBN(201 ,201) ,AYNP1(201) ,3YNP1(201) 
COMMON /INF2/ CCN ( 201 , 200 ) , CCT (201 , 200 ) , CYNPl ( 200 ) , CXNPl ( 200 ) 
COMMON /GUST/ UG(200 ) , VG(200 ) , XGF , UGUST , VGUST 
DIMENSION Bl(200) ,B2(200) ,AA(2) ,BB(2) 



RETRIEVE SOLUTION FROM A-MATRIX 

DO 50 I = l,NODTOT 

B1(I) = A(I,NP1) 

50 B2(I) = A(I,NP2) 

FIND VT AT TRAILING EDGE PANELS 
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DO 130 K 
IF (K .EQ, 
IF (K .EQ, 
XMID = 
YMID = 
VTANG 



AA(K) 

BB(K) 

DO 120 
AA(K) 
3B(K) 
CONTINUE 






= 1,2 

1 = 1 

I = NODTOT 
0.5 * (X(I) + X(I+1)) 

0.5 * (Y(I) + Y(I+1)) 

Ul. +UG ( I ) ) *COS ALF- VG ( I ) *S INALF+UX) *COSTHE 
+ ( ( 1 . +UG ( I ) ) *S INALF+VG ( I ) *COS ALF+UY ) *S INTHE 
+ OMEGA*(YMID*COSTHE(I) - XMID’^SINTHE(I) ) 

- AAN(I,NP1)*SS/WAKE 

VTANG + AAN(I,NP1)’^SS*GAMMA/WAKE 

= 1, NODTOT 

AA(K) + AAN(I,J) - BBN(I, J)*B1(J) 

3B(K) - 3BN(I, J)*B2(J) 




ADD CORE VORTEX CONTRIBUTION 



110 

100 

130 



IF (M .EQ. 1) GO TO 100 
MMl = M - 1 
DO 110 N = 1,MM1 
BB(K) = BB(K) + CCT(I,M)*CV(N) 
CONTINUE 
CONTINUE 
CONTINUE 



SATISFYING KUTTA CONDITION -- SOLVE FOR VORTEX STRENGTH 
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c 

c 

c 



Ec, 

FF 

GG 

RADI 

GAMK 




= SQRT(FF’^FF-EE*GG) 

= (-FF - RADI)/EE 

CALCULATE SOURCE STRENGTH 



- SS/DT 

+ 2.^SS’^GAim/DT 



DO 160 I = l.NODTOT 
160 QK(I) = GAMK^Bl(I) + B2(I) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

C SUBROUTINE NACA45(Z, THICK, CAMBER, BETA) C 

C C 

C EVALUATE THICKNESS AND CAMBER C 

C FOR NACA 4- OR 5-DIGIT AIRFOIL C 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE NACA45 ( Z , THI CK , CAMBER , BETA) 

COMMON /PAR/ NACA,TAU,EPSMAX,PTMAX 
THICK =0.0 

IF (Z .LT. l.E-10) GO TO 100 

THICK = 5.’^TAU^( .2969’^SQRT(Z) - Z’^(.126 + Z^(.3537 
+ - Z*(.2843 - Z^riOlS)))) 

100 IF (EPSMAX .EQ. 0.0) GO TO 130 
IF (NACA .GT. 9999) GO TO 140 
IF (Z .GT. PTMAX) GO TO 110 

CAMBER = EPSMAX/PTMAX/PTMAX^(2.*?TMAX - Z)^Z 
DCAMDX = 2. ’^EPSMAX/ PTMAX/ PTMAX’^( PTMAX - Z) 

GO TO 120 

110 CAMBER = EPSMAX/ (l.-PTMAX)’^’^2*(l. + Z - 2 . ^PTMAX) 1 . - Z) 

DCAMDX = 2. ’^EPSMAX/ (1. -PTMAX) ’^’^2* (PTMAX - 2) 

120 BETA = ATAN( DCAMDX) 

RETURN 

130 CAMBER =0.0 

BETA =0.0 

RETURN 

140 IF (Z .GT. PTMAX) GO TO 150 

W = Z/ PTMAX 

CAMBER = EPSMAX*W*((W - 3.)*W + 3. - PTMAX) 

DCAMDX = EPSMAX*3.^W*(1 . - W) /PTMAX 
GO TO 120 

150 CAMBER = EPSMAX*(1. - Z) 

DCAMDX = - EPSMAX 

GO TO 120 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

C SUBROUTINE PRESS (SINALF,COSALF, OMEGA, UX,UY) C 

c c 

C COMPUTE UNSTEADY FLOW PRESSURE DISTRIBUTION C 

C AND VELOCITY POTENTIAL AT MID-POINTS OF PANELS C 

C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE PRESS (5INALF , COSALF , OMEGA , UX , UY) 

COMMON /BOD/ ZFLAG , MLOWER , NUPP5R , NODTOT , X( 202 ) , Y ( 202 ) , 

+ COSTHE(201) ,SINTHE(201) ,SS,N?1,NP2 

COMMON /CPD/ CP (200) 

COMMON /NUM/ PI,PI2INV 

COMMON /SING/ 0(200) , GAMMA, QK( 200) , GAMK 
COMMON /WAK/ VYW,VXW,WAKE,DT 

COMMON /CORV/ CV( 200 K XC ( 200 ) , YC ( 200 ) ,M, TD , CCVX( 200 ) , CCVY( 200 ) 
COMMON /INFl/ AAN(201,201),BBN(201,201),AYNP1(201),BYNP1(201) 
COMMON /INF2/ CCN( 201 , 200 ) , CCT ( 201 , 200 ) , CYNPl ( 200 ) , CXNPl ( 200 ) 
COMMON /POT/ PHI(200) ,PHIK(200) 

COMMON /GUST/ UG(200 ) , VG( 200 ) , XGF , UGUST , VGUST 
COMMON /EXTV/ UE(200) 
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WRITE (6,1000) 

FIND TANGENTIAL VELOCITY VT AT MID-POINT OF I-TH PANEL 
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DO 130 

XMID 

YMID 

DX 

DY 

DIST 

VSX 

VSY 

VS- 

VTANG 



VTFRES 
VTANG 
DO 120 
VTANG 
CONTINUE 



OMEGA’^YMID + UX 
0MEGA*XMID + UY 



I = l,NODTOT 
= 0.5 (X(I) + X(I+1)) 

= 0,5 - Y(I) r Y(I+1)) 

= [y[t+iS yIiH 

= S0RT(DX’"DX+DY^DY) 

'l.+UG(I) )*COSALF-VG(I)’^SINALF 
' 1 . +UG ( I ) ) ’^SINALF+VGj( I ) ’^COSALF 
VSX’^VSX + VSY’^VSY 
= ( ( 1 . +UG ( I ) ) *C0SALF- VG ( I ) *S INALF+UX) *COSTHE (I) 

+ ( (^1 . +UG (.1 n *S INALF+VG ( I ) *COS ALF+UY ) INTHE ( I ) 
+ OMEGA’^(YMID’^COSTHE(I) - XMID^SINTHE ( I ) ) 

= VTANG 

= VTANG + SS’^(GAHMA-GAMK)*AAN(I,NP1)/WAKE 
J = 1,N0DT0T 

= VTANG - BBN(I, J)*0K(J) + AAN ( I , J ) *GAMK 



ADD CORE VORTEX CONTRIBUTION 

IF (M .EO, 1) GO TO 150 
HMl = M - 1 





DO 140 


N 


= 1,MM1 




VTANG 


— 


VTANG + CCT(I ,N)^CV(N) 


140 


CONTINUE 




150 


CONTINUE 








PHIK(I) 




(VTANG-VTFREE)*DIST 




C?(I) 




VS - VTANG’^VTANG 




UE(I) 




VTANG 


130 


CONTINUE 







COMPUTE DISTURBANCE POTENTIAL BY LINE INTEGRAL OF VELOCITY FIELD 

INTEGRATION FROM UPSTREAM (AT INFINITY) TO THE LEADING EDGE 

MPKI = 10 * NLOWER 

PINK =0.0 

XL =0.0 

DO 30 L = 1,NPHI 

FRACT = FLOAT (L) /FLOAT (NPHI) 

XLP = -10.0 * (1.0 - COS(0.5*PI*FRACT)) 

DELX = XL - XLP 

XMID = 0.5*(XL+XLP)*COSALF 

YMID = 0.5*(XL+XLP)*SINALF 

XL = XLP 

VELX = UGUST 

ADD CONTRIBUTION OF J-TH PANEL 

DO 20 J = 1,NP1 

DXJ = XMID - X(J) 

DXJP = XMID - X(J+1) 

DYJ = YMID - Y(J) 

DYJP = YMID - Y J+1) 

FLOG = .5^ALOG( (DXJP’^DXJP+DYJP=^DYJP)/(DXJ*DXJ+DYJ’^DYJ) ) 

FTAN = ATAN2(DYJP^DXJ-DXJP’^DYJ,DXJP’‘'DXJ+DYJP’^DYJ) 

CALMTJ = -COSALF*COSTHE(J) - SINALF’^SINTHE ( J) 

SALMTJ = -SINALF*C0STHE( J) + COSALF’^SINTHE ( J) 

APY = PI2INV*(FTAN*CALMTJ + FLOG’^SALMTJ ) 

BPY = PI2INV*(FL0G*CALMTJ - FTAN’^SALMT J) 

IF (J .EQ. NPl) GO TO 40 
VELX = VELX - BPY*QK(J) +GAMK’^APY 

GO TO 20 

40 VELX = VELX + SS’^APY’^ (GAMMA-GAMK) /WAKE 
20 CONTINUE 
C 
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ADD CORE VORTEX CONTRIBUTION 



I? (H .EQ. 1) GO TO 50 
= M - 1 
DO 60 N = 

DX = XMID - XC(N) 

DY = YMID - YC(N) 

DIST = SORT(DX’^DX-rDY^DY) 

COSTHN = DX/DIST 
SIHTHI'I = DY/DIST 

SALMTN = -SIMALF’^COSTHN + COSALF’^SINTHN 
CRT = -?I2II'IV^SALHTN/DIST 
60 VELX = VELX + CPT’^CV(N) 

50 CONTINUE 

PINK = PINK + VELX * DELX 
30 CONTINUE 

COMPUTE DISTURBANCE POTENTIAL AT MID-POINT OF I-TH PANEL 

LOWER SURFACE 

DO 230 I = l,NLOWER 
?H = -PINK 

DO 240 J = I,NL0WER 
240 ?H = ?H - ?HIX(J) 

?HIK(I) = ?H 
230 CONTINUE 

DO 270 I = l.MLOWER-l 

PHIKi'I) = 0.5^(PHIK(I) + ?HIK(I + 1)) 

270 CONTINUE 

?HIK(NLOWER) = 0 . 5^(PHIK(NLOWER) + PINK) 

UPPER SURFACE 



DO 250 I = NODTOT,NLOWER+1,-1 
PH = -PINK 

DO 260 J = NLOWER+1,1 
260 PK = ?H + PHIK(J) 

?HIK(I) = PH 
250 CONTINUE 

DO 280 I = NODTOT,NLOWER+2,-l 
PHIK(I) = 0.5^(PHIK(I) + PHIK(I-l)) 

280 CONTINUE 

PHIK(NLOWER+l) = 0 . 5*(PHIK(NL0WER+1 ) + PINK) 



COMPUTE CP AT MID POINT OF I-TH PANEL 



DO 290 I = 1,N0DT0T 

XMID = .5^(XU) + X(I+1)) 

CP(I) = CP(I) - 2,*(PHIK(I)-PHI(I))/DT 
WRITE (6,1050) I ,XMID,QK(I) ,GAMK,CP(I) ,UE(I) 

290 CONTINUE 

1000 F0RHAT(/,4X, ' J' ,4X, ' X( J) ' , 6X , ' Q( J) ' , 5X , 'GAMMA' ,5X, 
^ + 'CP(J)',6X,'V(J)',/) 

1050 FORMAT(I5,6F10.6) 

RETUPvN 



END 

CCCCCCCCCCCCCCCCCCCCCCCC 



:cccc( 






ccccccccc 






SUBROUTINE SETUP 



SETUP COORDINATES OF PANEL NODES AND SLOPES OF PANELS 
COORDINATES ARE READ FROM INPUT DATA FILE UNLESS 

THE AIRFOIL IS OF NACA XXXX OR NACA 230XX TYPE 



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE SETUP 

COMMON /BOD/ IFLAG , NLOWER ,NUPPER ,NODTOT , X( 202) , Y( 202 ) , 

+ COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON /MUM/ PI,PI2INV 
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onooonno 



onnnnn non non 



PI 

PI2INV 



= 3.1415926585 
= .5/PI 

SET COORDINATES OF NODES ON BODY SURFACE 



IF (IFLAG .NE. 0) GO TO 10 

NPOINT = NLOWER 

SIGN = -1.0 

NSTART = 0 

DO 110 NSURF =1,2 

DO 100 N = 1, NPOINT 

FRACT = FLOAT(N-l)/FLOAT(NPOINT) 

Z = .5^(1. - COS(PI’^FRACT)) 

I = NSTART + N 

CALL BODY(Z,SIGN,X(I) ,Y(I)) 

100 CONTINUE 

NPOINT = NUPPER 

SIGN =1.0 
NSTART = NLOWER 

110 CONTINUE 

NODTOT = NLOWER + NUPPER 
X(NODTOT+l) = X(l) 

Y(NODTOT+l) = Y(l) 

GO TO 20 

10 NODTOT = NLOWER + NUPPER 

READ (1,501) (X(I),I=l,NODTOT+l) 
WRITE (6,501) (X(I) ,I=l,NODTOT+l) 
READ (1,501) (Y(I),I=l,NODTOT+l) 
WRITE (6,501) (Y(I) ,I=l,NODTOT+l) 
501 FORMAT (6F10.6) 

20 NPl = NODTOT + 1 



200 



NP2 


= NODTOT + 2 


S 


ST SLOPES OF 


SS 


= 0.0 


DO 200 


I = 1, NODTOT 


DX 


= X(I+1) - X 


DY 


= Y(i+1) - Y 


DIST 


= SQRT(DX*DX 


SS 


= SS + DIST 


SINTHE(I 


) = DY/DIST 


costhe(i 

CONTINUE 

RETURN 


) = DX/DIST 



OF PANELS AND CALCULATE AIRFOIL PERIMETER 



END 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



SUBROUTINE TEWAK (SINALF , COSALF) 

COMPUTE WAKE ELEMENT AT THE TRAILING EDGE 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE TEWAK (SINALF , COSALF) 

COMMON /BOD/ IFLAG, NLOWER, NUPPER, NODTOT, X(202) ,Y(202) , 

+ COSTHE(201) ,SINTHE(201) ,SS,NP1,N?2 

COMMON /COF/ A(201 ,211) ,NEOS 
COMMON / S ING/ 0(200), GAMMA TOK ( 2 0 0 ) , GAMK 
COMMON /WAK/ VYW , VXW , W.AKE , DT 
COMMON /WAK2/ VYWK,VXWK 

COMMON /CORV/ CV(200 ) , XC ( 200 ) , YC ( 200 ) ,M , TD , CCVX(200) , CCVY(200 ) 
COMMON /IMFl/ AAN(201,201) ,BBN(201,201) ,AYNP1(201) ,BYNP1(201) 
COMMON /INF2/ CCN(201 , 200 ) , CCT ( 201 , 200 ) , CYNPl (200 ) , CXNPl (200 ) 
COMMON /GUST/ UG(200) ,VG(200) ,XGF,UGUST,VGUST 
XMID = 0.5 * (X(NPl) + X(NP2)) 

YMID = 0.5 * (y(np1) + y(nP2)) 

UGW =0.0 

VGW =0.0 

XG = XMID^COSALF + YMID’^SINALF 

IF (XG .GT. XGF) GO TO 10 
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nnnnn 



non non non nnnnnnnn nnn 



UGW = UGUST 

VGW = VGUST 

10 VYWK = (l.+UGW)’='SINALF+VGW’^COSALF 
VXWK = (l.+UGW)’^COSALF-VGW*SIMALF 

DO 120 J = 1,N0DT0T 

VYWK = VYWK + AYN?1(J)’^QK(J) + BYNPl ( J)’^GA2'IK 
120 VXWK = VXWK - BYNP1(J)^QK(J) + AYNP 1 ( J ) ’"GAMK 

ADD CORE VORTEX CONTRIBUTION 



IF (M .EQ. 1) GO TO 140 
MMl = M - 1 
DO 130 N = 1,MM1 
VYWK = VYWK + CYNP1(N)*CV(N) 
VXWK = VXWK + CXNP1(N)^CV(N) 
CONTINUE 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCi 
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140 



cccccccccccccccccccc 



SUBROUTINE VELDIS (SINALF , COSALF) 

COMPUTE STEADY FLOW PRESSURE DISTRIBUTION 

AND VELOCITY POTENTIAL AT MID-POINTS OF PANELS 



GCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE VELDIS (SINALF , COSALF ) 

COMMON /BOD/ IFLAG ,NLOWER ,NUPPER , NODTOT , X( 202 ) , Y( 202) , 

+ COSTHE(201) ,SINTHE(201) ,SS,NP1,NP2 

COMMON /COF/ A(201,211) ,KUTTA 
CCHMOM /CPD/ C?(200) 

COMMON /NUM/ PI,PI2INV 

COMMON /SING/ 0 (200 ), GAMMA , OK ( 200 ), GAMK 
COMMON /POT/ Pfil(200) ,PHIK(200) 

COMMON /GUST/ UG( 200 ) ,VG( 200 ), XGF , UGUST , VGUST 
COMMON /EXTV/ UE(200) 

WRITE (6,1000) 

RETRIEVE SOLUTION FROM A-MATRIX 



DO 50 I = 1, NODTOT 
50 Q(I) = A(I,KUTTA+1) 

GAMllA = A ( KUTT A , KUTT A+ 1 ) 

FIND VT AND CP AT MID-POINT OF I-TH PANEL 



DO 130 I = 1, NODTOT 

XMID = .5’^(X(I) + X(I+1)) 

YMID = .5*(Y(I) + Y(I+1)) 

VTANG = COSALF’^COSTHE(I) + SINALF’^SINTHE ( I ) 
VTFREE = VTANG 



ADD CONTRIBUTION OF J-TH PANEL 



100 



120 



DO 120 J = 1, NODTOT 
FLOG =0.0 

PTAN = PI 

:? J .EQ. :) GO TO 100 
DXJ - XMID - X(J) 

DXJP = XMID - X(J+1) 

DYJ = YMID - Y(J) 

DYJP = YMID - Y(J+1) 

FLOG = .5*AL0G({DXJP’^DXJP+DYJP’^DYJP)/(DXJ’^DXJ+DYJ*DYJ)) 

FTAN = ATAN2(DYJP’^DXJ-DXJP*DYJ,DXJP^DXJ+DYJP*DYJ) 

CTIMTJ = COSTHE(l)*COSTHE(J) + SINTHE ( I ) *SINTHE ( J) 

STIMTJ = SINTHE(I)*COSTHE(J) - COSTHE ( I ) ^SINTHE ( J) 

AA = PI2INV^(FTAN*CTIMTJ + FLOG^STIMTJ) 

B = PI2INV*(FLOG*CTIMTJ - FTAN’^STIMTJ) 

VTANG = VTANG - B*Q(J) +GAMMA’^AA 

CONTINUE 
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n o n n n n o 



non ooooo non ooooo ooo 



CP(I) = 1.0 - VTANG^VTANG 

UE ( I ) = VTANG 

WRITE (6,1050) ,GAMMA,CP(I) ,UE(I) 

INITIAL SET-UP FOR DISTURBANCE POTENTIAL CALCULATION 

DX = X(I+1) - X(I) 

DY = Y(I+1) - Y(I) 

DIST = SQRT(DX^DX+DY*DY) 

PHI (I) = (VTAI'IG-VTFREE)^DIST 

130 CONTINUE 



COMPUTE DISTURBANCE POTENTIAL BY LINE INTEGRAL OF VELOCITY FIELD 



INTEGRATION FROM UPSTREAM (AT INFINITY) TO THE LEADING EDGE 



NPHI 

PIN 

XL 

DO 30 

FRACT 

XLP 

DELX 

XMID 

YMID 

XL 

VELX 



= 10 ^ NLOWER 
= 0.0 
= 0.0 

L = 1,NPHI 

= FLOAT(L)/FLOAT(NPHI) 

= -10.0 (1.0 - COS(0.5*PI*FRACT)) 

= XL '- XLP 

= 0.5’^(XL+XLP)*COSALF 
= 0.5^(XL+XLP)’^SINALF 
= XLP 
= UGUST 



ADD CONTRIBUTION OF J-TH PANEL 



20 

30 



DO 20 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CALMTJ 

SALMTJ 

APY 

BPY 

VELX 

CONTINUE 

PIN 

CONTINUE 



= l,NODTOT 
XMID - X(J) 

XMID - X(J+1) 

YMID - Y(J) 

YMID - Y(J+1) 

. 5^ALOG( (DXJP’^DXJP+DYJP’^DYJP) / (DXJ*DXJ+DYJ^DYJ) ) 
ATAN2(DYJP*DXJ-DXJP’^DYJ,DXJP^DXJ+DYJP*DYJ) 
-COSALF’^COSTHE(J) - SINALF’^SINTHE ( J) 
-SINALF*COSTHE(J) + COSALF^SINTHE ( J) 
?I2INV’^(FTAN*CALMTJ + FLOG’^SALMTJ) 
?I2INV^(FLOG^CALMTJ - FTAN*SALMTJ) 

VELX - BPY’^Q(J) +GAMMA^APY 

PIN + VELX * DELX 



COMPUTE DISTURBANCE POTENTIAL AT MID-POINT OF I-TH PANEL 



LOWER SURFACE 

DO 230 I = 1, NLOWER 
PH = -PIN 

DO 240 J = I, NLOWER 
240 PH = PH - PHI(J) 

PHI ( I ) = PH 

230 CONTINUE 

DO 270 I = 1, NLOWER- 1 

?HI(I) = 0.5^(?HI(I) ^ PKI(I+1)) 

270 CONTINUE 

PHI (NLOWER) = 0.5^(PHI (NLOWER) + PIN) 
UPPER SURFACE 

DO 250 I = NODTOT,NLOWER+l,-l 
PH = -PIN 

DO 260 J = NLOWER+1,1 
260 PH = PH + PHI(J) 

PHI (I) = PH 

250 CONTINUE 

DO 280 I = NODTOT,NLOWER+2,-l 
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?HI(I) = 0.5^(?HI(I) r PHI(I-D) 

230 CONTINUE 

PHKNLOWER+l) = 0.5*(?HI(NL0WER+1) + PIN) 

1000 ?0P2-LaT(/,4X, ' J ' 4X, ' X( J) ' , 6X, ' Q ( J) ' , 5X, 'GAMMA' ,5X, 
+ 'C?(J)',5X,'V(J)',/) 

1050 FORMAT ( I 5, 5c iO. 6) 

RETURN 

END 
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APPENDIX B 

EXAMPLE INPUT DATA FOR PROGRAM U2DIIF 



THIS IS AN EXAMPLE OUTPUT DATA OBTAINABLE FROM PROGRAM U2DIIF 
AIRFOIL : MISES 8.4% THICKNESS (COORDINATES ARE INPUT BY USER) 

PANEL NUMBER : NLOWER = 25 , NUPPER =25 

AIRFOIL MOTION : MODIFIED RAMP AOA CHANGE ABOUT MID CHORD 

INITIAL AOA : 2.5 DEGREES 

FINAL AOA : 7.5 DEGREES 

AOA RISE TIME : 1.5 CHORD LENGTHS 

COMPUTATION TIME STEP ; 0.05 DURING TRANSIENT MOTION, INCREASES 

PROGRESSIVELY AFTER FINAL .^OA IS REACHED 



01 25 


25 










1.000000 


0.994858 


0.980866 


0.958884 


0.929536 


0.893455 


0.851308 


0.803815 


0.751753 


0.695948 


0.637271 


0,576620 


0.514918 


0.453098 


0.392084 


0.332794 


0.276105 


0,222365 


0.173361 


0.129819 


0.091393 


0.059146 


0.033560 


0.015010 


0.003767 


0.000000 


0.003767 


0.015008 


0.033560 


0.059146 


0.091393 


0.129819 


0.173861 


0.222865 


0.276105 


0,332791 


0.392082 


0.453095 


0.514915 


0.576617 


0.637266 


0.695946 


0.751750 


0.803815 


0.851308 


0.893455 


0.929536 


0.958884 


0.980866 


0.994858 


1.000000 








0.000000 


-0.000732 


-0.002784 


-0.005721 


-0.009351 


-0.013459 


-0.017837 


-0.022235 


-0.026618 


-0.030671 


-0.034289 


-0.037341 


-0.039712 


-0.041314 


-0.042083 


-0.041979 


-0.040979 


-0.039096 


-0.036360 


-0.032820 


-0.028555 


-0.023651 


-0.018220 


-0.012379 


-0.006259 


0.000000 


0.006259 


0.012379 


0.018220 


0.023651 


0.028555 


0.032820 


0.036360 


0.039096 


0.040979 


0.041979 


0.042083 


0.041314 


0.039712 


0.037341 


0.034289 


0.030671 


0.026618 


0.022285 


0.017837 


0.013459 


0.009351 


0.005721 


0.002784 


0.000782 


0.000000 








2.50000 


5.000000 


1.5 


o 

o 


0.5 


O 

o 


0.000000 


0.000000 


0.000000 








2.000000 


0.050000 


0.0001 


0.000000 
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APPENDIX C 

EXAMPLE OUTPUT DATA FROM PROGRAM U2DIIF 



DATA READ FROM FILE CODE 1 



11 



THIS IS AN EXAMPLE OUTPUT DATA OBTAINABLE FROM PROGRAM U2DIIF 
AIRFOIL : MISES 8.4% THICKNESS (COORDINATES ARE INPUT BY USER) 
PANEL NUMBER : NLOWER = 25 , NUPPER = 25 

AIRFOIL MOTION : MODIFIED RAMP AO A CHANGE ABOUT MID CHORD 
INITIAL AOA : 2.5 DEGREES 
FINAL AOA : 7.5 DEGREES 



AOA RISE TIME ; 1.5 CHORD LENGTHS 

COMPUTATION TIME STEP : 0.05 DURING TRANSIENT MOTION, INCREASES 

PROGRESSIVELY AFTER FINAL AOA IS REACHED 



1 25 

1.000000 
0.851308 
0.514913 
0.173361 
0.003767 
0.091393 
0.392082 
0.751750 
0.980866 
0.000000 
-0.017837 
-0.039712 
-0.036360 
-0.006259 
0.028555 
0.042083 
0.026618 
0.002784 
2.500000 
0.000000 

2.000000 



25 

0.994858 

0.803815 

0.453098 

0.129819 

0.000000 

0.129319 

0.453095 

0.803815 

0.994858 

-0.000782 

-0.022285 

-0.041314 

-0.032820 

0.000000 

0.032820 

0.041314 

0.022285 

0.000782 

5.000000 

0.000000 

0.050000 



0.980866 

0.751753 

0.392084 

0.091393 

0.003767 

0.173861 

0.514915 

0.851308 

1.000000 

-0.002784 

-0.026618 

-0.042083 

-0.028555 

0.006259 

0.036360 

0.039712 

0.017837 

0.000000 

1.500000 

0.000000 

0.000100 



0.958884 

0.695948 

0.332794 

0.059146 

0.015008 

0.222865 

0.576617 

0.893455 

-0.005721 

-0.030671 

-0.041979 

-0.023651 

0.012379 

0.039096 

0.037341 

0.013459 

0.000000 

0.000000 



0.929536 

0.637271 

0.276105 

0.033560 

0.033560 

0.276105 

0.637266 

0.929536 

-0.009351 

-0.034289 

-0.040979 

-0.018220 

0.018220 

0.040979 

0.034289 

0.009351 

0.500000 



0.893455 

0.578620 

0.222865 

0.015010 

0.059146 

0.332791 

0.695946 

0.958884 

-0.013459 

-0.037341 

-0.039096 

-0.012379 

0.023651 

0.041979 

0.030671 

0.005721 

0.000000 



0.000000 



AIRFOIL PERIMETER LENGTH = 2.018599 



STEADY FLOW SOLUTION AT ALPHA = 2.500000 





x(J) 


Q(J) 


G.AMMA 


C?(J) 


V(J) 


1 


0.997429 


0.355723 


0.074003 


0 . 316305 


-0.326859 


2 


0.987862 


0.356105 


0.074003 


0.206074 


-0.391025 


3 


0.969875 


0.365026 


0.074003 


0.133790 


-0.930704 


4 


0.944210 


0.378836 


0.074003 


0.082276 


-0.957979 


5 


0.911495 


0.394973 


0.074003 


0.043033 


-0.978247 


6 


0.872381 


0.412926 


0.074003 


0.012034 


-0.993965 


7 


0.827561 


0.432568 


0.074003 


-0.012724 


-1.006342 


8 


0.777784 


0.453710 


0.074003 


-0.032414 


-1.016078 


9 


0.723850 


0.476112 


0.074003 


-0.048010 


-1.023724 


10 


0.666609 


0.500047 


0.074003 


-0.059905 


-1.029517 


11 


0.606945 


0.525455 


0.074003 


-0.068405 


-1.033637 


12 


0.545769 


0.552654 


0.074003 


-0.073563 


-1.036129 
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13 


0.484008 


0.581715 


0.074003 


-0.075309 


-1.036971 


14 


0.422591 


0.612382 


0.074003 


-0.073531 


-1.036114 


15 


0.362439 


0.546480 


0.074003 


-0.067928 


-1.033406 


16 


0.304449 


0.683297 


0.074003 


-0.057668 


-1.028430 


17 


0.249485 


0.723595 


0.074003 


-0.041891 


-1.020731 


18 


0.198363 


0.763523 


0.074003 


-0.019114 


-1.009512 


19 


0.151840 


0.319732 


0.074003 


0.013374 


-0.993290 


20 


0.110606 


0.379132 


0.074003 


0.059669 


-0.969707 


21 


0.075269 


0.951277 


0.074003 


0.127774 


-0.933931 


22 


0.046353 


1.043426 


0.. 074003 


0.232984 


-0.875795 


23 


0.024285 


1.172784 


0.074003 


0.409236 


-0.768612 


24 


0.009388 


1.380641 


0.074003 


0.727179 


-0.522323 


25 


0.001884 


1.653644 


0.074003 


0.945112 


0.234282 


26 


0.001884 


0.545367 


0.074003 


-0.815326 


1.347341 


27 


0.009387 


-0.275210 


0.074003 


-1.084184 


1.443670 


28 


0.024284 


-0.497235 


0.074003 


-0.872601 


1.368430 


29 


0.046353 


-0.580394 


0.074003 


-0.723499 


1.312821 


30 


0.075269 


-0.613032 


0.074003 


-0.624167 


1.274428 


31 


0.110606 


-0.636699 


0.074003 


-0.552954 


1.246176 


32 


0.151840 


-0.645594 


0.074003 


-0.498371 


1.224080 


33 


0.198363 


-0.649755 


0.074003 


-0.453794 


1.205734 


34 


0.249485 


-0.650971 


0.074003 


-0.415592 


1.189787 


35 


0.304448 


-0.650596 


0.074003 


-0.381557 


1.175397 


36 


0.362436 


-0.649624 


0.074003 


-0.349957 


1.161877 


37 


0.422538 


-0.648020 


0.074003 


-0.319875 


1.148858 


38 


0.484005 


-0.646281 


0.074003 


-0.290579 


1.136037 


39 


0.545766 


-0.644572 


0.074003 


-0.261352 


1.123099 


40 


0.606941 


-0.642990 


0.074003 


-0.231604 


1.109776 


41 


0 . 666606 


-0.641396 


0.074003 


-0.200941 


1.095875 


42 


0.723848 


-0.639980 


0 .074003 


-0.168763 


1.081094 


43 


0.777782 


-0.538504 


0.074003 


-0.134640 


1.065195 


d4 


0 . 327 561 


-0.537214 


0.074003 


-0.097676 


1.047701 


45 


0.372381 


-0.535330 


0.074003 


-0.056864 


1.028039 


46 


0.911495 


-0.534250 


0.074003 


-0.010900 


1.005436 


47 


0.944210 


-0.632196 


0.074003 


0.042413 


0.978564 


48 


0.969875 


-0.629793 


0.074003 


0.107609 


0.944665 


49 


0.987862 


-0.625960 


0.074003 


0.193062 


0.898297 


50 


0.997429 


-0.619834 


0.074003 


0.316307 


0.826858 



CD = 0.000829 CL = 0.303076 CM = -0.080325 



'k'k'k'k'k:k'k:k:k:k'k:k'k:k:k:k'k'k'k'k:k'k:k'k'k'k:k'k:k'k'k'k:k:k'k:k:k:k'k 

*** BEGIN UNSTEADY FLOW SOLUTION **** 

'k'k'k'k'k:k'k:k:k7^:k:k:k:k7Z7^:k:ki^'k:ki^i^'k:k:ky^:k7^i^:k:k:k:k:k:k:k:k'k 



TIME STEP TK = 0.050000 



TK - TKMl = 0.050000 



ALPHA(T) = 2.516295 

U(T) = 0.000000 



OMEGA(T) = -0.011248 

V(T) = -0.005624 



NITR VXW 



'7YW WAKE THETA GAMMA 



0 0.999048 0.043619 

1 0.907832 0.005991 

2 0.904138 0.007297 

3 0.903985 0.007241 



0.050000 0.043633 
0.045393 0.006600 
0.045208 0.008070 
0.045201 0.008010 



0.740032E-01 

0.744799E-01 

0.744662E-01 

0.744652E-01 



CONVERGED SOLUTION OBTAINED AFTER NITR = 3 

J X(J) Q(J) GAMMA CP(J) V(J) 

1 0.997429 0.435837 0.074466 0.299470 -0.838202 
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2 


0.987362 


0.422352 


0.074466 


0.196120 


-0.899569 


3 


0.969375 


0.421793 


0.074466 


0.132077 


-0.936976 




0.944210 


0.427256 


0.074466 


0.088439 


-0.962416 


5 


0.911495 


0.435974 


0.074466 


0.055859 


-0.981181 


G 


0.372381 


0.447181 


0.074466 


0.029322 


-0.995676 


7 


0.327561 


0.460598 


0.074466 


0.008149 


-1.007037 


8 


0.777784 


0.475913 


0.074466 


-0.010442 


-1.015962 


9 


0.723850 


0.492830 


0.074466 


-0.026773 


-1.022942 


10 


0 • 666609 


0.511571 


0.074466 


-0.041132 


-1.028238 




0,606945 


0.532044 


0.074466 


-0.053494 


-1.032001 


12 


0.545769 


0.554564 


0.074466 


-0.063531 


-1.034262 


13 


0.484008 


0.579201 


0.074466 


-0.070990 


-1.035018 


14 


0.422591 


0.606206 


0.074466 


-0.075185 


-1.034204 


15 


0.362439 


0.635915 


0.074466 


-0.075487 


-1.031672 


16 


0.304449 


0.669133 


0.074466 


-0.070722 


-1.027008 


1 n 


0.249485 


0.706140 


0.074466 


-0.059745 


-1.019773 


13 


0.198363 


0.748100 


0.074466 


-0.040826 


-1.009171 


19 


0.151340 


0.796633 


0.074466 


-0.011154 


-0.993762 


20 


0.110606 


0.853325 


0.074466 


0.033396 


-0.971231 


21 


0.075269 


0.924099 


0.074466 


0.100715 


-0.936853 


. 22 


0.046353 


1.014820 


0.074466 


0.205876 


-0.880676 


23 


0.024235 


1.143358 


0.074466 


0.382654 


-0.776542 


24 


0.009338 


1.351849 


0.074466 


0.703606 


-0.535871 


25 


0.001884 


1.634712 


0.074466 


0.952740 


0.209596 


40 


0.001834 


0.564272 


0.074466 


-0.746972 


1.322641 


27 


0.009387 


-0.246433 


0.074466 


-1.037466 


1.430099 


28 


0.024234 


-0 . 467318 


0.074466 


-0.838061 


1.360476 


29 


0.046353 


-0.551795 


0.074466 


-0.693563 


1.307915 


30 


0.075269 


-0.590860 


0.074466 


-0.596473 


1.271481 


31 


0.110606 


-0.611392 


0.074466 


-0.527117 


1.244626 


32 


0.151340 


-0.622549 


0.074466 


-0.474318 


1.223580 


33 


0,193363 


-0.529333 


0.074466 


-0.433321 


1.206047 


34 


0.249435 


-0 . 633515 


0.074466 


-0.399064 


1.190716 


35 


0.304448 


-0 , 636433 


0.074466 


-0.369826 


1.176793 


36 


0.362436 


-0.639059 


0.074466 


-0.343641 


1.163586 


37 


0.422588 


-0.641343 


0.074466 


-0.319346 


1.150745 


38 


0.484005 


-0.643768 


0.074466 


-0.295850 


1.137963 


39 


0.545766 


-0.646483 


0.074466 


-0.272088 


1.124936 


40 


0.606941 


-0.649573 


0.074466 


-0.247068 


1.111386 


41 


0 . 666606 


-0.652913 


0.074466 


-0.220048 


1.097122 


42 


0.723848 


-0.656699 


0.074466 


-0.190134 


1.081849 


43 


0.777782 


-0.660707 


0.074466 


-0.156552 


1.065285 


44 


0.827561 


-0.665236 


0.074466 


-0.118313 


1.046980 


45 


0.872381 


-0.670135 


0.074466 


-0.074279 


1.026307 


46 


0.911495 


-0.675261 


0.074466 


-0.023233 


1.002476 


47 


0.944210 


-0.680614 


0.074466 


0.036802 


0.974108 


43 


0.969875 


-0.686543 


0.074466 


0.109850 


0.938386 


49 


0.987862 


-0.692719 


0.074466 


0.203356 


0.889792 


50 


0.997429 


-0.699982 


0.074466 


0.333160 


0.815602 


CD 


= 0.001539 CL = 


0.302054 CM = 


-0.088450 


TRAILING VORTICES DATA 








M 


X(M) 


Y(M) 


CIRC 






1 


1.022599 


0.000131 


-0.000933 







TIME STEP TK = 0.749999 



TK - TKMl = 0.050000 



ALPHA(T) = 4.999996 

U(T) = 0.000000 



OMEGA(T) = -0.087266 

V(T) = -0.043633 



NITR VXW VYW 

0 0.905684 -0.000916 



WAKE THETA GAMMA 

0.045284 -0.001012 0.103235E+00 
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1 0.905735 -0.000649 0.045287 -0.000717 



0.106563E+00 



CONVERGED SOLUTION OBTAINED AFTER NITR = 1 



J 


X(J) 


Q(J) 


GAMMA 


CP(J) 


V(J) 


1 


0.997429 


1.115997 


0.106565 


0.311649 


-0.908159 


2 


0.987862 


1.060333 


0.106565 


0.221106 


-0.957033 


3 


0.969875 


1.031653 


0.106565 


0.170306 


-0.983709 


4 


0.944210 


1.013085 


0.106565 


0.141163 


-0.998976 


5 


0.911495 


0.998141 


0.106565 


0.123949 


-1.008098 


6 


0.872381 


0.985265 


0.106565 


0.114073 


-1.013416 


7 


0.827561 


0.974241 


0.106565 


0.109016 


-1.016142 


8 


0.777784 


0.964993 


0.106565 


0.107136 


-1.017013 


9 


0.723850 


0.957451 


0.106565 


0.107130 


-1.016594 


10 


0.666609 


0.952120 


0.106565 


0.108458 


-1.015089 


11 


0.606945 


0.949235 


0.106565 


0.110657 


-1.012669 


12 


0.545769 


0.949393 


0.106565 


0.113665 


-1.009317 


13 


0.484008 


0.952960 


0.106565 


0.117540 


-1.004971 


14 


0.422591 


0.960451 


0.106565 


0.122480 


-0.999505 


15 


0.362439 


0.972455 


0.106565 


0.128967 


-0.992654 


16 


0.304449 


0.990059 


0.106565 


0.138078 


-0.983855 


17 


0.249485 


1.013807 


0.106565 


0.150962 


-0.972486 


18 


0.198363 


1.045081 


0.106565 


0.169616 


-0.957451 


19 


0.151840 


1.085785 


0.106565 


0.197364 


-0.936848 


20 


0.110606 


1.133034 


0.106565 


0.239070 


-0.907675 


21 


0.075269 


1.206453 


0.106565 


0.303822 


-0.863894 


22 


0.046353 


1.298127 


0.106565 


0.407818 


-0.793054 


23 


0.024285 


1.428849 


0.106565 


0.583190 


-0.663132 


24 


0.009388 


1.631804 


0.106565 


0.872165 


-0.369686 


25 


0.001884 


1.819807 


0.106565 


0.765241 


0.485826 


26 


0.001884 


0.373179 


0.106565 


-1.559307 


1.595828 


27 


0.009387 


-0.529374 


0.106565 


-1.564167 


1.590937 


28 


0.024284 


-0.755145 


0.106565 


-1.202155 


1.468068 


29 


0.046353 


-0.836350 


0.106565 


-0.991312 


1.389589 


30 


0.075269 


-0.874103 


0.106565 


-0.864650 


1.338450 


31 


0.110606 


-0.896239 


0.106565 


-0.781035 


1.302194 


32 


0.151840 


-0.912096 


0.106565 


-0.721021 


1.274529 


33 


0.198363 


-0.926607 


0.106565 


-0.674005 


1.251843 


34 


0.249485 


-0.941343 


0.106565 


-0.634257 


1.232132 


35 


0.304448 


-0.957406 


0.106565 


-0.598321 


1.214148 


36 


0.362436 


-0.975544 


0.106565 


-0.563539 


1.196886 


37 


0.422588 


-0.995441 


0.106565 


-0.528592 


1.179831 


38 


0.484005 


-1.017302 


0.106565 


-0.492391 


1.162518 


39 


0.545766 


-1.041010 


0.106565 


-0.454043 


1.144540 


40 


0.606941 


-1.066387 


0.106565 


-0.412924 


1.125555 


41 


0 . 666606 


-1.093023 


0.106565 


-0.368740 


1.105335 


42 


0.723848 


-1.120818 


0.106565 


-0.321026 


1.083543 


43 


0.777782 


-1.149241 


0.106565 


-0.269645 


1.059939 


44 


0.827561 


-1.178310 


0.106565 


-0.213999 


1.034027 


45 


0.872381 


-1.207644 


0.106565 


-0.153563 


1.005269 


46 


0.911495 


-1.236895 


0.106565 


-0.087684 


0.972965 


47 


0.944210 


-1.265976 


0.106565 


-0.014749 


0.935775 


48 


0.969875 


-1.296105 


0.106565 


0.069103 


0.890819 


49 


0.987862 


-1.330154 


0.106565 


0.171240 


0.832327 


50 


0.997429 


-1.380203 


0.106565 


0.308161 


0.746084 


CD 


= 0.033887 CL = 


: 0.64533i 


8 CM - 


-0.224298 


TRAILING VORTICES DATA 








M 


X(M) 


Y(M) 


CIRC 






1 


1.700760 


0.077052 


-0.000933 






2 


1.651289 


0.072156 


-0.001625 






3 


1.602002 


0.066331 


-0.002229 






4 


1.552845 


0.060031 


-0.002784 






5 


1.503779 


0.053471 


-0.003304 






6 


1.454797 


0.046788 


-0.003797 






7 


1.405895 


0.040107 


-0.004256 
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8 


1.357057 


0.033554 -0.004687 


9 


1.308305 


0.027216 -0.005100 


10 


1.259669 


0.021159 -0.005481 


11 


1.211178 


0.015541 -0.005801 


12 


1.162916 


0.010481 -0.006100 


13 


1.115048 


0.006079 -0.006349 


14 


1.067918 


0.002406 -0.006559 


15 


1.022643 


-0.000016 -0.006720 



TIME STEP TK = 


1.449992 




TK - TKMl 


= 0.050000 


AL?KA(T) = 7. 


483698 


OMEGA (T) 


= -0.011 


249 


U(T) 


= 0. 


000000 


V(T) 


= -0.005625 


NITR 


VXW 


VYW 


WAKE 


THETA 


GAMMA 


0 


0.901699 


0.009872 


0.045088 


0.010948 


0.145377E+00 


1 


0.901200 


0.011028 


0.045063 


0.012236 


0.146997E+00 


COMVERGED SOLUTION OBTAINED AFTER NITR = 1 




J 


X(J) 


Q(J) 


GAMMA 


CP(J) 


V(J) 


1 


0.997429 


1.101898 


0.146996 


0.332408 


-0.864078 


2 


0.987862 


1.099609 


0.146996 


0.228004 


-0.921138 


3 


0.969875 


1.116333 


0.146996 


0.160643 


-0.955077 


4 


0.944210 


1.140985 


0.146996 


0.114863 


-0.976579 


0 


0.911495 


1.163619 


0.145996 


0.082627 


-0.990894 


6 


0.872381 


1.197904 


0.146996 


0.060523 


-1.000257 


7 


0.827561 


1.228671 


0.146996- 


0.046710 


-1.005378 


8 


0.777784 


1.260826 


0.146996 


0.039936 


-1.008490 


9 


0.723850 


1.294203 


0.146996 


0.039112 


-1 . 008660 


10 


0.666609 


1.329222 


0.146996 


0.043766 


-1.006560 


11 


0.606945 


1.366027 


0.146996 


0.053332 


-1.002352 


12 


0.545769 


1.405159 


0.146996 


0.067653 


-0.995947 


13 


0.484008 


1.446882 


0.146996 


0.086578 


-0.987214 


14 


0.422591 


1.491619 


0.146996 


0.110103 


-0.975912 


15 


0.362439 


1.539856 


0.146996 


0.138557 


-0.961606 


16 


0.304449 


1.592619 


0.146996 


0.172964 


-0.943442 


17 


0.249485 


1.650390 


0.146996 


0.214399 


-0.920465 


18 


0.198363 


1.714437 


0.146996 


0.265094 


-0.890942 


19 


0.151840 


1.786607 


0.146996 


0.328746 


-0.851951 


20 


0.110606 


1.868882 


0.146996 


0.410545 


-0.798852 


21 


0.075269 


1.965495 


0.146996 


0.519597 


-0.722328 


22 


0.046353 


2.082446 


0.146996 


0.668455 


-0.603478 


23 


0.024285 


2.230950 


0.146996 


0.866926 


-0.394782 


24 


0.009388 


2.419835 


0.146996 


1.009466 


0.050859 


25 


0.001884 


2.339782 


0.146996 


-0.471799 


1.214887 


26 


0.001884 


-0.156346 


0.146996 


-4.389847 


2.320095 


27 


0.009387 


-1.322127 


0.146996 


-3.033496 


2.003043 


28 


0.024284 


-1.560176 


0.146996 


-2.015200 


1.727204 


29 


0.046353 


-1.622657 


0.146996 


-1.505832 


1.569752 


30 


0.075269 


-1.634560 


0.146996 


-1.212790 


1 .4705^9 


31 


0.110606 


-1.623102 


0.1^6996 


-1.021304 


1.401554 


32 


0.151340 


-1.613625 


0.146996 


-0.385616 


1.350007 


33 


0.193363 


-1.596423 


0.146996 


-0 . 730655 


1.309004 


34 


0.249485 


-1.578180 


0.146996 


-0.695048 


1.274895 


35 


0.304448 


-1.560042 


0.146996 


-0.621833 


1.245407 


36 


0.362436 


-1.542861 


0.146996 


-0.556431 


1.218916 


37 


0.422588 


-1.526391 


0.146996 


-0.496616 


1.194569 


38 


0.484005 


-1.510876 


0.146996 


-0.440592 


1.171595 


39 


0.545766 


-1.496317 


0.146996 


-0.387199 


1.149437 


40 


0.606941 


-1.482638 


0.146996 


-0.335611 


1.127617 


41 


0.666606 


-1.469497 


0.146996 


-0.285518 


1.105863 


42 


0.723848 


-1.456874 


0.146996 


-0.236273 


1.083745 


43 


0.777782 


-1.444335 


0.146996 


-0.187542 


1.060991 
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44 


0.327561 


-1.431955 


0.146996 


-0.138453 


1.037087 


45 


0.372381 


-1.419472 


0.146996 


-0.038061 


1.011468 


46 


0.911495 


-1.406547 


0.146996 


-0.035055 


0.983411 


47 


0.944210 


-1.393024 


0.146996 


0.023027 


0.951546 


48 


0.969875 


-1.379863 


0.146996 


0.091342 


0.912885 


49 


0.987862 


-1.368515 


0.146996 


0.178655 


0.861777 


50 


0.997429 


-1.365097 


0.146996 


0.303827 


0.784388 


CD 


= 0.0309 


'56 CL ^ 


= 0.713821 


II 

u 


-0.190685 


TRAILING VORTICES DATA 








M 


X(M)- 


Y(M) 


CIRC 






1 


2.383780 


0.232008 


-0.000933 




■ 


2 


2.333846 


0.225080 


-0.001625 






3 


2.284404 


0.216393 


-0.002229 






4 


2.235273 


0.206695 


-0.002784 








2.136347 


0.196301 


-0.003304 






6 


2.137600 


0.185376 


-0.003797 






1 


2.089004 


0.174067 


-0.004256 






8 


2.040442 


0.162550 


-0.004687 






9 


1.991957 


0.150853 


-0.005100 






10 


1.943572 


0.138989 


-0.005481 






11 


1.395155 


0.127193 


-0.005301 






12 


1.846653 


0.115534 


-0.006100 






13 


1.798126 


0.104170 


-0.006349 






14 


1.749496 


0.093078 


-0.006559 






15 


1.700793 


0.082359 


-0.006720 






16 


1.651988 


0.072087 


-0.006839 






17 


1.603077 


0.062348 


-0.006896 






18 


1.554079 


0.053196 


-0.006916 






19 


1.505019 


0 .044648 


-0.006885 






20 


1.455917 


0.036753 


-0.006795 






21 


1.406791 


0.029591 


-0.006648 






22 


1.357685 


0.023172 


-0.006443 






23 


1.308651 


0.017529 


-0.006177 






24 


1.259734 


0.012684 


-0.005852 






25 


1.211021 


0.008642 


-0.005465 






26 


1.162620 


0.005401 


-0.005014 






27 


1.114706 


0.002948 


-0.004498 






28 


1.067624 


0.001210 


-0.003917 






29 


1.022530 


0.000276 


-0.003269 







TIME STEP TK = 1.999990 



TK - TKMl = 0.200000 



ALPHA(T) = 7.500000 

U(T) = 0.000000 



OMEGA(T) = 0.000000 

V(T) = 0.000000 



NITR VXW 



VYW WAKE 



THETA GAMMA 



0 0.939783 

1 0.948367 

2 0.948647 



0.026143 0.138029 
0.030228 0.139770 
0.030081 0.189825 



0.027811 0.156187E+00 
0.031363 0.160749E+00 
0.031699 0.160765E+00 



CONVERGED SOLUTION OBTAINED AFTER NITR = 2 



J X(J) 



Q(J) GAMMA CP(J) V(J) 



1 0.997429 

2 0.987862 

3 0.969875 

4 0.944210 

5 0.911495 

6 0.872381 

7 0.827561 



1.116135 0.160766 
1.110748 0.160766 
1.126034 0.160766 
1.150497 0.160766 
1.179177 0.160766 
1.210700 0.160766 
1.244766 0.160766 



0.320772 -0.851401 
0.221351 -0.907726 
0.158555 -0.941336 
0.116479 -0.962935 
0.086724 -0.977640 
0.065709 -0.987588 
0.051593 -0.993866 
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3 


0.777784 


1.281123 


0.160766 


0.043212 


-0.997146 


9 


0.723850 


1.319423 


0 . 160765 


0.039694 


-0.997912 


10 


0 . 566609 


1.359943 


0.160766 


0.040813 


-0.996296 


11 


0.506945 


1.402723 


0.160766 


0.046325 


-0.992422 


' 7 


0 . 545769 


1.448176 


0.160766 


0.056472 


-0.986157 


13 


0.484008 


1.496434 


0.150756 


0.071440 


-0.977368 


14 


0.422591 


1.547997 


0.160766 


0.091673 


-0.965767 


15 


0.362439 


1.503142 


0.160766 


0.117858 


-0.950894 


^ -- 

_o 


0.304449 


1.562397 


0.160766 


0.151404 


-0.931843 


17 


0.249435 


1.727721 


0.160756 


0.193674 


-0.907603 


IS 


0.198353 


1.793850 


0.160766 


0.247145 


-0.876338 


19 


0.151840 


1.378113 


0.160756 


0.315652 


-0.834967 


20 


0.110606 


1.967481 


0.160766 


0.404369 


-0.778575 


21 


0.075269 


2.071099 


0.160766 


0.522025 


-0.697329 


22 


0.046353 


2.194821 


0.160766 


0.679633 


-0.571292 


23 


0.024285 


2.349161 


0.160756 


0.381013 


-0.350465 


24 


0.009383 


2.539111 


0.160756 


0.937480 


0.119060 


25 


0.001384 


2.420488 


0.160766 


-0.772858 


1.331881 


26 


0.001384 


-0.236994 


0.160766 


-4.940613 


2.437122 


27 


0.009387 


-1.441366 


0.160766 


-3.294956 


2.071290 


23 


0.024284 


-1.678374 


0.160766 


-2.145366 


1.771569 


29 


0.046353 


-1.735029 


0.160766 


-1.575494 


1.601988 


30 


0.075269 


-1.740151 


0.160766 


-1.248166 


1 .495590 


31 


0.110606 


-1.726694 


0.150766 


-1.035354 


1.421868 


32 


0.151340 


-1.705142 


0.160766 


-0.884674 


1.367020 


33 


0.198363 


-1.580851 


0 . 160766 


-0.770207 


1.323621 


■54. 


0.249435 


-1.655530 


0.150766 


-0.678837 


1.287744 


35 


0.304443 


-1.530341 


0.160766 


-0.602871 


1.256978 


36 


0.362435 


-1.606176 


0.160766 


-0.537036 


1.229567 


37 


0.422538 


-1.582302 


0.160766 


-0.478581 


1.204600 


38 


0.434005 


-1.560513 


0.160765 


-0.425234 


1.131278 


39 


0.545765 


-1.539373 


0.160766 


-0.375544 


1.158993 


40 


0.506941 


-1.519334 


0.160766 


-0.327563 


1.137220 




0 . 566605 


-1.500277 


0 . 160766 


-0.281143 


1 . 1 15668 


42 


0.723848 


-1.482141 


0.160766 


-0.235015 


1.093871 


43 


0.777782 


-1.464664 


0.160766 


-0.188536 


1.071527 


44 


0.827561 


-1.448072 


0.160766 


-0.140512 


1.048044 


45 


0.372381 


-1.432243 


0.160766 


-0.089803 


1.022803 


46 


0.911495 


-1.417014 


0.160766 


-0.035089 


0.995022 


47 


0.944210 


-1.402392 


0.160766 


0.026050 


0.963238 


43 


0.969875 


-1.389331 


0.160766 


0.098463 


0.924430 


49 


0.987362 


-1.379196 


0.160766 


0.190332 


0.873000 


50 


0.997429 


-1.379105 


0.160766 


0.319847 


0.795184 


CD 


= 0.022037 CL = 


= 0.709635 CM = 


-0.185327 


AILING VORTICES DATA 








M 


X(M) 


Y(M) 


CIRC 






1 


2.923984 


0.318366 


-0,000933 






2 


2.873369 


0.311788 


-0.001625 






3 


2.823534 


0.302817 


-0.002229 






4 


2.774192 


0.292439 


-0.002784 






5 


2.725182 


0.281064 


-0.003304 






6 


2.576492 


0.263830 


-0.003797 






7 


2.528036 


0.255054 


-0.004256 






3 


2.579755 


0.342831 


-0.004637 






9 


2.531597 


0.229197 


-0.005100 






10 


2.433743 


0.215114 


-0.005431 






11 


2.435836 


0.200928 


-0.005801 






12 


2.387918 


0.186791 


-0.006100 






13 


2.339993 


0.172653 


-0.006349 






14 


2.291922 


0.158710 


-0.006559 






15 


2.243772 


0.144992 


-0.006720 






16 


2.195477 


0.131596 


-0.006839 






17 


2.146968 


0.118650 


-0.006896 






13 


2.098253 


0.106217 


-0.006916 






19 


2.049394 


0.094288 


-0.006885 






20 


2.000362 


0.082951 


-0.006795 
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21 


1.951119 


0.072314 


22 


1.901708 


0.062412 


23 


1.852145 


0.053298 


24 


1.802427 


0.045045 


25 


1.752599 


0.037677 


26 


1.702674 


0.031252 


27 


1.652678 


0.025824 


28 


1.602603 


0.021481 


29 


1.552431 


0.018390 


30 


1.501920 


0.017193 


31 


1.450588 


0.016109 


32 


1.378677 


0.012411 


33 


1.258439 


0.007189 


34 


1.094864 


0.003008 



I 



- 0.006648 

- 0.006443 

- 0.006177 

- 0.005852 

- 0.005465 

- 0.005014 

- 0.004498 

- 0.003917 

- 0.003269 

- 0.002553 

- 0.002761 

- 0.005504 

- 0.007734 

- 0.009244 
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